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FOREWORD 


This report constitutes final documentation of work performed by per- 
sonnel of Lockheed's Huntsville Research & Engineering Center for NASA- 
Langley Research Center on Contract NAS1-15783. The purpose of the report 
is the presentation of the GIM code modifications and the latest flowfield 
calculations performed with the General Interpolants Method (GIM) computer 
code. The work describes the GIM code reacting flow model, the hyperbolic 
steady-state Euler version of the code, investigation of a solution adaptive 
grid algorithm, and improvements to the turbulence models. The authors 
gratefully acknowledge the contribution to and continued support of this 
work by the NASA-Langley Contract Monitors J. L. Hunt and R.C. Rogers. 
Inquiries concerning this report should be directed to: 


John F. Stalnaker 

Lockheed-Huntsville Research & Engineering Center 
4800 Bradford Drive 
Huntsville, Alabama 35807 

Telephone: (205) 837-1800, ext. 401 
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1. INTRODUCTION AND SUMMARY 


The General Interpolants Method (GIM) code was developed to analyze co- 
mplex flow fields which defy solution by simple methods. The code uses num- 
erical difference techniques to solve the full three-dimensional time- 
averaged Navier-Stokes equations in arbitrary geometric domains. The nu- 
merical analogs of the differential equations are derived by representing 
each flow variable with general interpolant functions. The point of de- 
parture then requires that a weighted integral of interpolants be zero over 
the flow domain. By choosing the weight functions to be the interpolants 
themselves, the GIM formulation can produce the classical implicit differ- 
ence schemes. Choosing the weight functions to be orthogonal to the inter- 
polant functions produces explicit finite difference type discrete analogs. 
By appropriate choice of constants in the weight functions, the GIM becomes 
analogous to standard finite difference schemes such as centered, backward, 
forward, windward and multi-step predictor-corrector schemes. The GIM 
analogs, however, are automatically produced for arbitrary geometric flow 
domains and hence are a general point of departure and provide flexibility 
in the choice of differencing schemes. 

The GIM computer code was originally written for the CDC 7600 machine. 
The first effort that was accomplished under Contract to NASA-Langley was 
the conversion and reprogramming of the code for the CDC-STAR (now termed 
the CYBER 200 series) vector processor. The GIM-STAR code was then exer- 
cised for three-dimensional exhaust flows for application to Scramjet engine 
studies. The next sequential study in this computational fluid dynamics 
effort consisted of the development and application of a parabolized GIM 
algorithm, computation of the flow, including spillage, in a model aircraft 
inlet and investigation of linearized block implicit schemes for GIM appli- 
cation. These tasks were accomplished under two contracts through a 
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cooperative effort of the Hypersonic Aerodynamics and Hypersonic Propulsion 
Branches. These efforts then continued to include program modification to 
fully utilize the features of the CYBER 203, inclusion of the implicit 
MacCormack algorithm, inclusion of algebraic and differential equation tur- 
bulence models, creation of an interactive input program, calculation of the 
inviscid flow about a wing-body configuration, and the incorporation of a 
finite rate nonequilibrium reacting flow model. 

The most current effort, which is the subject of this report, is a 
continuation of the GIM code development and application on the CYBER 203 
machine. Objectives of this effort include the following. 

• Complete the development of the finite rate chemically reacting 
flow model. 

• Validate the reacting flow model for hydrogen-air and 
hydrocarbon-air combustion by analyzing an experimental case. 

• Develop a hyperbolic marching Euler solver for inviscid flow- 
field analysis. 

• Investigate solution-adaptive grid algorithms for inclusion in 
the GIM code. 

• Complete the formulation of the turbulence models in the GIM 
code. 

Certain of these stated objectives were accomplished in full and other 
partially. The plan of attack was to organize a set of overall tasks, some 
of which are interrelated, aimed at meeting the major objectives. This 
report is organized into sections with each section independently presenting 
details of the major tasks. The following is a summary of these sections: 

THE GIM CODE CHEMICALLY REACTING FLOW MODEL 


The GIM code chemically reacting flow model has been developed to 
facilitate the calculation of a broad spectrum of flow fields involving 
chemical reactions. The GIM code model allows the calculation of frozen 
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flow, equilibrium flow, and finite-rate nonequilibrium flow. Both an ex- 
plicit and an implicit time integration scheme are available at the option 
of the user. All of the chemistry calculations are tailored to be readily 
compatible with the overall G1M code methodology and are highly vectorized 
so as to obtain maximum computational efficiency on the CDC CYBER 200 series 
of supercomputers. Furthermore, the model is completely general in that any 
number of reactions of the general form 


N 



i=l 
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£ 



A. 
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j - 1 
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Number of 
Reactions 


can be considered. In addition to the regular GIM code integration input, 
the user must input the equation describing each reaction, reaction rate 
data, and tables of the thermodynamic properties of each species as a func- 
tion of temperature. An initial estimate of the species mass fraction dis- 
tribution is also required. As output, the GIM code produces the regular 
integration output plus species mass fraction distributions. 

THE GIM HYPERBOLIC STEADY EULER SOLVER 

The General Interpolants Method solution methodology was brought to 
bear on a spatial-marching MacCormack algorithm for the solution of the 
hyperbolic steady Euler equations of inviscid flow. The resulting code has 
the computational efficiency resulting from the combination of the explicit 
MacCormack algorithm, the progressive assembly of the GIM difference analogs 
(see Appendix D), and the CYBER 200 series of supercomputers. Further 
features of the method include the geometric versatility of all GIM deriv- 
atives and dynamic monitoring of the marching step size with the ability to 
interpolate secondary cross planes in order to maintain stability. 
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INVESTIGATION OF A SOLUTION ADAPTIVE GRID ALGORITHM FOR THE GIM CODE 


A study was made to determine the feasibility and suitability of a 
solution adaptive grid scheme for the GIM code. This task consisted of a 
literature survey of currently available methods and development of a one- 
dimensional code implementing a test algorithm. The approach taken was to 
be the "boundary conforming" coordinate transformation whereby all the grid 
movement is contained in the form of metric coefficients. The grid veloc- 
ities were determined by setting them proportional to the local pressure 
gradient. The one-dimensional algorithm and an example result are presented 
in this report. The example case is a standing shock wave in a tube. The 
grid points are shown to migrate toward the discontinuity in pressure. The 
shock profile is seen to improve over the fixed grid case in that the 
"wiggles" due to MacCormack differencing are no longer present. We conclude 
that adaptive grid methods are feasible for inclusion in the GIM code and 
that a three-dimensional study should be started toward this goal. 

TURBULENCE MODELING IN THE GIM CODE 

Section 5 describes improvements made to the turbulence models in the 
GIM code. The Baldwin- Lomax algebraic eddy viscosity model was extended to 
include two-dimensional and axisymmetric regions bounded by more than one 
wall. It was further modified to allow the user to limit the region of 
application of the model. The Turbulent Kinetic Energy (TKE) two-equation 
model was completed with the inclusion of wall function boundary conditions. 

RECOMMENDATIONS 


An important part of any research effort is the publishing of the 
findings and results. In addition to NASA Contractor Reports, the results 
of the GIM code development and applications have been presented at many 
professional society meetings and published in the open literature. The 
authors would like to acknowledge the realization of this aspect of research 
by the Langley personnel and thank them for their support in these efforts 
to disseminate knowledge. 
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Another equally important aspect of computational research is the 
effort to keep current and potential users of computer codes abreast of 
changes and modifications to the code. To this end a GIM Code User's Bul- 
letin has been circulated as warranted and three courses have been given to 
aquaint users with the latest versions of the code. 

To further enhance the GIM code's capability and to utilize its current 
potential, the following actions are recommended. 

• Include the new streamlined assembly procedure in all GIM INTEG 
modules . 

• Investigate solution-adaptive methods including adaptive grids, 
time steps, and difference schemes. 

• Include higher order interpolants into the present geometry pr- 
ogram and include other state-of-the-art grid generators. 

• Restructure the code for storage optimization. 

• Investigate matching techniques to allow the code to dynamic- 
ally switch among the hyperbolic, parabolic, and elliptic 
options. 

• Use the GIM approach as a solution methodology to take advan- 
tage of its geometric and algorithmic consistency for the 
development of better solution techniques. 

Finally, the use of the GIM code by Langley personnel in their many 
research efforts is encouraged and recommended. Only through use and users' 
input can a code with the capability of the GIM code be developed along 
rational lines. 
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2. GIM CODE CHEMICALLY REACTING FLOW MODEL 


2 . 1 INTRODUCTION 


The GIM Code chemically reacting flow model has been developed to 
facilitate the calculation of a broad spectrum of flow fields involving 
chemical reactions. The GIM code model allows the calculation of frozen 
flow, equilibrium flow, and finite rate nonequilibrium flow. Both an ex- 
plicit and an implicit time integration scheme are available at the option 
of the user. All of the chemistry calculations are tailored to be readily 
compatible with the overall GIM code methodology and are highly vectorized 
so that maximum computational efficiency can be obtained on the CDC CYBER 
200 series of supercomputers. Furthermore, the model is completely general 
in that any number of reactions of the general form 
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can be considered. In addition to the regular GIM code Integration input, 
the user must input the equation describing each reaction, reaction rate 
data, and tables of the thermodynamic properties of each species as a func- 
tion of temperature. An initial estimate of the species mass fraction dis- 
tribution is also required. As output, the GIM code produces the regular 
Integration output plus species mass fraction distributions. 

This section details the theory and algorithms employed in the GIM code 
chemically reacting flow model and describes the input required to implement 
the model. Section 2.2 presents the governing equations used in the model, 
while Section 2.3 describes the computational model. Section 2.4 details 
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the frozen flow and equilibrium flow models incorporated in the general 
chemically reacting flow model. Section 2.5 is an input guide describing 
use of the model as well as the inputs required. A sample nonequilibrium 
finite rate reacting chemistry calculation is presented in Section 2.6. 
Nomenclature and references are presented in Sections 2.7 and 2.8, respec- 
tively. 

2.2 GOVERNING EQUATIONS 


The governing equations for general three-dimensional chemically react- 
ing flow with N species and M reactions as modeled in the GIM code include 
the following: 

Global Continuity Equation 
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Momentum Equations 
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Energy Equation 
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Species Continuity Equations 
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Equation of State 
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Written in the conservation vector format employed in the GIM code, Eqs, 
(2.1) through (2.4) become 
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where U, Ej, and H are given below: 
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2.3 COMPUTATIONAL MODEL 


2.3.1 General 

The integration of Eqs. (2.1) through (2.4), or equivalently equation 
(2.6), is not as straightforward as it is when the flow field is considered 
alone. The species continuity Eqs. (2.4a) contain an algebraic production 
or source term, namely w^ which most often involves a characteristic time 
scale much smaller than the characteristic time scale typically associated 
with the flowfield integration. This results in a set of "stiff" partial 
differential equations which can be computationally difficult to integrate. 
Two techniques are available in the GIM code to integrate this set of equa- 
tions. Both techniques employ the MacCormack (Ref. 1) explicit predictor- 
corrector scheme to integrate the flow field Eqs* (2.1) through (2.3). The 
species continuity Eqs. (2.4) are integrated using an iterative predictor- 
corrector-chemistry scheme in which the chemistry terms can be solved either 
explicitly or implicitly at the option of the user. The explicit scheme can 
be subject to severe time step restrictions but can be computationally quite 
efficient. The implicit scheme is not nearly so restricted with regard to 
the time step but requires more computational effort per time step. 

2.3.2 Flowfield Differencing 

The MacCormack explicit predictor-corrector differencing scheme is 
always used to integrate the flowfield Eqs. (2.1) through (2.3), or altern- 
atively the first 5 of Eqs. (2.6). This scheme as used in the GIM code is 
given by: 
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- z £ <*r> 

u 2 

(2.8a) 

B[6U n+1 ] 

-r dU^ 1 

(2.8b) 

u n+1 = 

u ° + T (6ljn+1 + <su n+1 ) 

(2.8c) 

at 

U n +! + bt ({U n+1 - <5U n+1 ) 



2.3.3 Chemistry Differencing 


The species continuity Eqs. (2.4) or alternately the last N of Eqs. 
(2.6) are integrated using a predictor-corrector-chemistry scheme in which 
the chemistry production terms can be solved either explicitly or im- 
plicitly, at user option. This method employs the explicit MacCormack 
predictor-corrector differencing scheme to predict the influence of convec- 
tion and diffusion which is then used in the final iterative chemistry step. 
This scheme as used in the GIM code is given below: 
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where $ is an adjustable parameter which can vary between zero and one, 

0 < * < 1. The value of 4> determines to what degree the chemistry produc- 
tion terms are treated explicitly or implicitly. A value of <J> « 0 renders 
the chemistry terms fully explicit, while a value of = 1 yields a fully 
implicit treatment. Typically, a value of 0 = 0.5 is used which yields a 
second order trapezoidal-like integration scheme. If the magnitudes of 
the chemistry production terms are so large that the species continuity 
equations become stiff, then 4> is adjusted toward a value of 1.0 so that the 
equations are integrated in a nearly fully implicit mode. 

Equations (2.11) are implicit in U n+1 through the H n+1 term. These 
equations can be solved in an explicit manner through linearization of the 

and the method of successive substitutions or they can be solved as 
is via the Newton-Raphson method for nonlinear equations. 

2. 3. 3.1 Explicit Solution Technique 

Equations (2.11) can be solved explicitly by linearizing the 
term and then using the method of successive substitutions. This is done as 
follows : 
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where (3H/3U) is the Jacobian matrix associated with the H vector and where 
superscript (£) indicates evaluation using the £ iteration approximation 
of U n+ ^. Combine Eqs. (2.11) and (2.12) to obtain 
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At 
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At 
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Equation (2.13) represents N equations which are solved explicitly at each 
node in the flow field for the quantity 
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The values of 
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Equations (2.13) and (2.14) are iterated on £ beginning with 
- U n until the difference becomes negligible, at which 

time U n+1 = U (£+1) . 


2. 3. 3. 2 Implicit Solution Technique 


Equations (2.11) can be solved implicitly by using the Newton-Raphson 
method for nonlinear equations as follows: 
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Combine Eqs. (2.11) and define F 
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where superscript (£) indicates evaluation using the iteration approx- 
imation of U n+ "*' and 
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Combining Eqs. (2.15) yields 



+ (1 - (j)) H n + <j>H (£) 


(2.16) 


where I is the identity matrix and (3H/3U) is the Jacobian matrix associated 
with the H vector. Equation (2.16) is a system of N equations which must be 
solved at each node in the flow field for the quantity 

r u a+1 > - U ( *\ 
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The values of are given by 

U <t + D , „<*> + 4t , ° a+1) A - U K) ) (2.17, 

Equations (2.16) and (2.17) are iterated on H beginning with ^ 

= U n until the difference becomes negligible, at which 

time U n+1 - U^ +1 \ 

2.3.4 Decode Procedures 

After obtaining the conserved variables at the new time level from Eqs. 
(2.7) and (2.8) plus Eqs. (2.9) through (2.11), the primitive variables p, 
u^, U£, u^, £, c^, p, and T must be decoded. From the definition of 
U given in Eq. (2.6b) it can be seen that most of the primitive variables 
can be obtained simply by dividing U by p. However, a more complex rela- 
tionship exists between p, T, £, and c^ so that a special decode procedure 
must be devised in order to obtain p and T. This decode procedure as ap- 
plied in the GIM code is given as follows 

Given: p, u^ , U2, u^, £> and c^ i = N after each step 

Find: p.T 
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(2.18c) 


'flow field (T) " ^~i q2 + p = ^--| q2 + RT 


and solve implicitly for T and then p = p RT. 


To solve Eq. (2.18) directly, define h ± (T) 
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With the value of the function h^(T) available at J discrete temperature 
abscissa, i.e., 


W " h i,j for Tj = (j - 1) AT j = 1,J 

the value of the function at any temperature point in the range T 1 < T < 

Tj can be represented by a second order Lagrangian interpolation poly- 
nomial, i.e., 
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where j is selected such that 


i T - Tj| < AT j - 2 J-l 

This can be further simplified using the constant AT assumption: 
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Equation (2.20) is a simple quadratic equation in (T/AT) of the form 
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with the solution given by 
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2.3,5 Vectorizatlon Procedures 


As stated earlier, all of the GIM code chemistry procedures are highly 
vectorized so as to obtain maximum computational efficiency on the CYBER 200 
series of supercomputers. This is true for the production and rate Eqs. 
(2.4b) through (2.4g); the explicit and implicit integration Eqs. (2.7) 
through (2.11), (2.12) through (2.14), and (2.15) through (2.17); and the 
decode procedure Eqs. (2.18) through (2.22). The solution of the implicit 
chemistry Eqs. (2.16) which involves N simultaneous equations at each node 
in the flow field is also vectorized. This is possible because of the alge- 
braic form of the production term H with the result that Eq. (2.16) contains 
unknown values only at the node of interest and does not involve terms at 
any other node in the flow field. Therefore, if there are NN nodes in the 
flow field under consideration, there will be NN systems of N simultaneous 
equations to solve each system of which is independent of all other systems 


• 21e) 
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but identical in form. Since the same sequence of mathematical operations 
must he applied to each of the NN systems of equations in order to solve 
them, this sequence of operations can be applied once to vector descriptors 
pointing to corresponding entries in the NN coefficient matrices and right- 
hand side, as opposed to applying the same sequence to each of the NN sys- 
tems. These alternatives are illustrated below: 

Let Eq. (2.16) be represented as 


where 


m - i-mQ 


M K = RHS 


V = 


,U+D „(*) 


i £<■;>• t A 

L j-i 3 j=i 3 


(e" + 1 ) + <1— 4>) H n + <}>H W 


Vectorized Approach 


L5 










Scalar Approach 



RHS 1 


M mT 

NN 

X 

NN 


rhs nn 


2.3.6 Artificial Numerical Diffusion 


Limited experience has indicated that some form of artificial numerical 
diffusion or viscosity must be introduced into the species continuity equa- 
tions (2.4a) to keep the solution smooth and physically consistent. This is 
particularly so in region of large concentration gradients on course compu- 
tational grids. The form for this artificial numerical diffusion term was 
developed following the motivation of McRae, Goodin, and Seinfeld (Ref. 2-2) 
and Forester (Ref. 2-3). The form consists of a simple second order 
diffusion-like term with a variable coefficient that responds to species 
gradients: 


ts [- 

. 1-1 2 


I Ac . | 

l 


(p u j ta j> si: 


(2.24) 


where 


Ac, 


(P 4k.) 


a global coefficient of 0(1); 

9cj 

“ i ii 

is a normalized gradient-sensitive 


Ax . - — , so that 
2 9 x j 


Acj 

=7 


coefficient; and 

= coefficient for consistency with convection terms. 


This artificial numerical diffusion is incorporated into the species conti- 
nuity Eqs . (2.4a) through the diffusion coefficient 2)y 
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where 







j 


real diffusion coefficient 


artificial diffusion 
i Ac, 


u . Ax 
3 j 


coefficient 


2.4 FROZEN FLOW AND EQUILIBRIUM FLOW MODELS 
2.4.1 Reaction Rate Limits 


Frozen flow and equilibrium flow are by definition the flow produced in 
the limit as the reaction rates become respectively infinitesimally small 
and infinitely large. Since the reaction rates appear explicitly only in 
the species continuity Eqs. (2.4a) via the production term the frozen 
flow and equilibrium flow limits can be obtained by considering the effect 
on these equations as the reaction rates become infinitesimally small and 
infinitely large. 


Substituting Eqs. (2.4e) into Eq. (2.4b) yields 


w. 


M 

E 


A v y k 


v! * v" 

N /pc„\ £j x 

- — n 


n 


i~l\ m. 


C j * =1 ' m Jt 


(2.25a) 


M 


. . a k f, x j 

3=1 3 
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where 


X 


j 



(2.25b) 


Then substituting Eq. (2.25) into the species continuity Eqs. (2.6) 
yields 




(2.26) 


from which the frozen flow and equilibrium flow limits can be determined. 
2.4.2 Frozen Flow 


In the frozen flow limit as the reaction rates become infinitesimally 
small the effect of chemical production becomes negligible so that convec- 
tion and diffusion dominate the time rate of change of the species concen- 
trations. In this limit the species continuity Eqs. (2.26) become 



0 (Frozen Flow) 


(2.27) 


Since w = 0 and hence H = 0, the chemistry step of the differencing scheme, 
Eq. (2.11), is superfluous and is not used. 


2.4.3 Equilibrium Flow 


The equilibrium flow limit can be obtained by dividing Eq. (2.26) by each 
, one at a time, and allowing to become infinitely large. This limiting 

j j 

process applied for the M different k f ’ s yields 
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0 


j I > 2 , . . . , 


(2,28a) 


X. 


J 


M 


from which 

w = H = 0 (2.28b) 

which is the steady state equilibrium flow limit* 


In order to obtain this steady state limit, the unsteady terra in Eq. 
(2.26) is retained. Furthermore, for stability considerations, the arti- 
ficial numerical diffusion term is also retained. Thus the equilibrium flow 
limit of the species continuity equations as used in the GIM code is given 
by Eqs. (2.6) with the convection terms omitted, i.e., 


3t 


(p c ± ) + 



- 0 


i=l,N 


(2.29) 


(Equilibrium Flow) 


This limiting form of the species continuity equations allows the equilib- 
rium solution to be obtained wherein the chemical production dominates the 
time rate of change of the species concentrations with convection and real 
diffusion being negligible. 

2.5 INPUT GUIDE 

2.5.1 General 

The GIM code chemically reacting flow model consists of a set of UPDATE 
directives and code to be used with the INTEG integration module and the 
DYNDIM dynamic dimensioning module. These UPDATE directives are located on 
a semi-private file named CHEMOD on user number 838700C. CHEMOD contains 
two logical records. The first record contains the UPDATE directives for 
the INTEG integration module; the second contains the UPDATE directives for 
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the DYNDIM dynamic dimensioning module. Both sets must be used with their 
respecitve modules in order to generate the GIM code chemically reacting 
flow model. 

2.5.2 NOS Side Runstream Information 


Figure 2-1 illustrates a typical NOS "f ront-side" or "Z-machine" run- 
stream for use with the GIM code chemically reacting flow model. 

2.5.3 CYBER 203 Runstream Information 


Figure 2-2 illustrates a typical CYBER 203 runstream for use with the 
GIM code chemically reacting flow model. 

2.5.4 Dynamic Dimensioning Input Data Summary 

The GIM code chemically reacting flow model requires dynamic dimension- 
ing input data that differ somewhat from that normally associated with 
INTEG. This section presents a summary of these data and a brief descrip- 
tion of each input variable. 


The dynamic dimensioning data for the GIM code chemically reacting flow 
model are input on one card in a 915 integer format. The integers are 
always right justified in five colume increments. The input variables are: 

MN , I DIM , NSPEC , NREAC t NC ATS , NSP , MNB , METHOD , IDYN 
MN 

The total number of nodes in the problem (elliptic run) or the number 
of nodes in the initial cross plane (quasi-parabolic run). 
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/JOB 

/NOSEQ 

RUNCHM,CM60000,T400. 

USER( 7 5097 8C , PASWRD ) 

CHARGE (01 2 34 56 ,LRC) 
GET(OLDPL=INTEG > CHEMOD/UN=838700C) 

UPDATE/ F,C=TAPE8 , L=A12 ) 

RETURN ( OLDP L ) 

GET(OLDPL=UTILOPL/UN=838700C) 

UPDATE (I=CHEMOD) 

FTN ( I “COMPILE , B=DYNDIM ,L=0 , A) 

MAP (OFF) 

DYNDIM. 

RENAME ( INTGS=TAPE3 ) 

RETURN ( OLDPL, DYNDIM .CHEMOD , TAPE8 , COMPILE ) 
ATTACH(FILE20) 

ATTACH(FILE17 or FILE21) 

TOSTAR(INPUT,INTGS,FILE20,FILE17=BI or FILE21=BI) 
DAYFILE ( NOSDAY) 

REPLACE (NOSDAY) 

EXIT. 

DA YFILE( NOSDAY) 

REPLACE (NOSDAY) 

EXIT. 

/EOR 

*READ, CHEMOD 


Any other UPDATE Directives 


/EOR 

6534 2 5 2 0 329 326 2 0 

/EOR 


Fig. 2-1 - NOS Runstream 
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STORE 750978 400SDS RUNCHM B 

STRSIDE ,T500 . 

REQUEST (FILE2 2/ 2000) 

REQUEST ( JANAF/2000 ) 

FORTRAN (I =INTGS , B=INTGB/ 300 , L=ERRFTN, OPT=B ) 

LOAD (INTGB , OU=ERRLOD , CN=INTEGO ,4000 , CDF =5000 

,GRLP=*TMVEC , *PRIM, *EBUF, *VPR0P ,GRLP=*UBUF, *BOUND , *DELXYZ , *AXSYM 
,GRLP=*XBUF1 , GRLP=*XBUF2 ,GRLP=*TAUP , *TAU, *TAUF, *QPNOD, *SECORD 
,GRSP=*IOUNIT, *CNTRL, *TRANSP , *TDATA, *VECP , *SQ , *PM, *SUBSBC, *USER 
, MZONE , *QPCOM, *QPPRNT, *CVGCOM, *CHEM, *THERMO 
,GR0L=*Q3MAP) 

INTEGO. 

T0S(Z=7 50978C > FILE22, JANAF) 

DAYFILE(STRDAY) 

EXIT. 

TOAS (Z=7 50978C,ERRFTN,ERRLOD,FILE22 , JANAF) 

DAYFILE(STRDAY) 

EXIT. 

/EOR 


INTEG Input Data 


Chemically reacting flow model input data 


/EOR 


Fig. 2-2 - CYBER 203 Runstream 
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I DIM 


The spatial dimension of the problem. 


IDIM = 1 
= 2 
= 3 


axisymmetric flow 
two-dimensional flow 
three-dimensional flow 


NSPEC 


The number of reactive chemical species in the problem. 


NREAC 


The number of chemical reaction mechanisms (equations) considered. 


NCATS 


The number of nonreacting catalytic or "third-body" species. 


NSP 


The number of special node terms. 


MNB 


The number of boundary node terms. 


METHOD 


Elliptic/Quasi-Parabolic flag. 


METHOD £ 2 
> 2 


Elliptic run 
Quasi-Parabolic run 
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IDYN 


Dynamic assembly flag. 

IDYN = 0 Regular GIM code run 

^ 0 Dynamic assembly GIM code run 

2.5.5 Inpu t Ca rd Summa ry 

The GIM code chemically reacting flow model requires some chemistry 
input data in addition to and immediately following the normal INTEG inte- 
gration input data. This section presents a summary of the input cards and 
formats used to input this chemistry data. A description of each input 
variable and its options is presented in Section 2.5.6. After a user be- 
comes familiar with the chemistry inputs, this summarized input guide can be 
used to quickly identify each card and its contents. 

Several formats are used to input the chemistry data to the GIM code 
chemically reacting flow model. These include: 

ALPHANUMERIC A1 , A6, and A8 

INTEGER II and 15 

REAL F5.0 and E10.0 

Integers are always right justified in one or five column increments. 

Decimal or floating point data occupy five or ten columns each with, prefer- 
ably, a decimal point punched on the card. 

The GIM code chemically reacting flow model is activated by selecting 
the appropriate value for the INTEG input variable ISPEC (see Ref. 2-4) as 
follows : 
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ISPEC = 0 regular GIM code single gas calculation 
(no chemistry) 

= 1 frozen flow calculation 

= 2 equilibrium flow calculation 

= 3 nonequilibrium finite rate reaction calculation 
For ISPEC > 0, the input data shown on the following page must be provided: 
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CARD TYPE 


VARIABLE LIST/FORMAT 


Cl 


C2 


NSPEC, NRSPC, NREAC, NCATS, NTIP, ITCHEM, ISTIFF, IXOPT, 
I OPT, IMOLE, ICHUN 

(1115) 


CSTOL(I) , 1-1,5 

(5E10.0) 


C3 


C4 


C5 


C6 


C7 


C8 


C9 


CIO 


Cll 


CPCF, HCF, SCF, RRCF (LCONV(I), 1-1,4) 
(4E210.0.4A8) 


LSPID(I) , WM(I), HFM(I) 

(A6,4X,2E10.0) 


[TIP(J) , (CPHS(J ,K,I) ,K=1,3) , TIP( J+l) , 
(CPHS(J+1,K, I) ,K=1,3] J=1,NTIP,2 
(8E10.0) 


LSPID(I+NSPEC) 

(A6) 

WFH(I, J) , J-l, NSPEC 

(16F5.0) 


1=1, NCATS 


1 


1=1, NSPEC 


t(RTYPE(I,J) ,1=1,2) , (LREAC(I),I=1, 78) J J-l, NREAC 
(211, 78A1) 


[(RRATC(I,K, J) ,1=1,5] K-1,1 or 2 J=l, NREAC 
(5E10.0) 

NJ, INC, NTOT, ITYPE 
(415) 

► 

CS(I), 1=1, NSPEC 

(8E10.0) 


[ ] denotes card is repeated one or more times as indicated. 

{ } denotes card group is repeated one or more times as indicated. 
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2.5.6 Description of Input Data 


This section presents a description of the input variables listed in 
Section 2.5.5. Each variable is identified as to its usage in the GIM code 
chemially reacting flow model with options and standard values also given. 
All 11 card types are not necessarily input for a given problem. Some of 
the control variables on card type Cl, for example, dictate which options 
have been selected and hence which input cards are required. This informa- 
tion is given in the discussion of each variable to be input. Each input 
card which is read is also printed out to aid the user in debugging a 
problem setup. 

CARD TYPE Cl - Chemistry Control Variables 

FORMAT (815) 


NSPEC 

The total number of independent chemical species in the problem being 

run. 

NSPEC £ 1 

Any number of species may be considered in a given problem, subject only to 
the limitations of storage and run time (CRUs) available. A value of NSPEC 
less than 1 is an error and terminates the run immediately. 

NRSPC 

The number of reactive chemical species (as opposed to inert species) 
in the problem being run. 

1 < NRSPC £ NSPEC 

NRSPC should in general be set equal to NSPEC. However, if inert chemical 
species are present, then setting NRSPC less than NSPEC can save substantial 
amounts of computational resources. 
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NREAC 


The number of chemical reaction mechanisms (equations) considered. 

NREAC £ 1 chemically reacting flow 

£ 0 frozen flow (no reactions) 

Any number of reactions may be considered in a given problem, again subject 
only to the limitations of storage and run time. 

NCATS 


The number of nonreacting catalytic or "third-body" species. 

NCATS £ 0 

Any number of catalytic species may be included in a given problem. 


NTIP 


The number of thermodynamic interpolation points in the species thermo- 
dynamic data tables. The user must input a set of thermodynamic data tables 
for each of the NSPEC chemical species in the problem. These tables consist 
of the specific heat at constant pressure, sensible enthalpy, and entropy at 
each of NTIP temperature points. These data can either be read in as input 
or can be read from a previously created coded file named JANAF, depending 
on the sign of NTIP. 


NTIP > 0 read data from input 

NTIP < 0 read data from JANAF file 

If NTIP is greater than 0, the JANAF file Is automatically created for use 
on the next run. A value of NTIP equal to 0 is an error and terminates the 
run immediately. 
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ITCHEM 


This variable controls frequency of the calls to the chemistry rou- 
tines. The chemistry routines are called every ITCHEM c iteration (time 
step). Generally, for frozen flow (ISPEC = 1) or nonequilibrium flow (ISPEC 
= 3) ITCHEM should be set to 1 so that the chemistry routines are called on 
every iteration. For equilibrium flow (ISPEC = 2) or for flows where the 
chemistry is near steady state ITCHEM can be set to larger values. A value 
of ITCHEM equal to 0 deactivates the chemistry routines so that they are 
never called. In this case the species distributions will remain at their 
initial values. 

ISTIFF 

Stiff Chemistry Equation Option. If the magnitude of the species pro- 
duction tern, w (and hence H) , becomes large, a numerical solution of the 
species continuity equations can be difficult to obtain for all but the 
smallest of time steps. In order to overcome this severe time step restric- 
tion, the parameter <f) can be adjusted toward a value of 1.0. This renders 
the treatment of the chemistry terms more implicit and easier to solve for 
larger time steps at the expense of losing second order accurte tracking of 
the transient behavior. Control of the parameter <J> is selected via ISTIFF. 
i 

i 

j ISTIFF = 0 nonstiff equations; cj> = 0.5 for second order 

' accurate trapezoidal-like integration of 

chemistry production terms. 

= 1 stiff equations; 0.5 <_ $ <_ 1.0, calculated 

internally so that chemistry production terms 
are treated more implicitly as |w| increases. 
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IXOPT 


Implicit/explicit chemistry integration option. The user can select 
the implicit chemistry solution technique, the explicit chemistry solution 
technique or a combination of the two through the use of IXOPT: 

IXOPT “ 1 explicit chemistry solution technique 

= 0 implicit/explicit chemistry solution technique 
= -1 implicit chemistry solution technique 

If the user selects the combined explicit/implicit chemistry solution tech- 
nique (IXOPT - 0), the program determines which technique to use on each 
time step based on the maximum local change of the species concentrations 
and the CSTOL tolerances (card type C2). Frozen flow (ISPEC == 1) automatic- 
ally selects IXOPT = 1 for the explicit technique only. 

I OPT 


Chemical species output option. The user can select the type of chem- 
ical species data to be output with the usual flowfield data via I0PT: 


I0PT = 0 flow field only; no chemistry data 

* 1 flow field + species mass (or mole) 

fractions 

- 2 flow field + species mass (or mole) 
fractions + species production rates 


IMOLE 


i 


i 


Mole fraction/mass fraction option. The GIM code chemically reacting 
flow model inputs and outputs the species concentrations in terms of mass 
fraction (by default). Input and output in terms of mole fractions can be 
selected using IMOLE as indicated below: 
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IMOLE = 0 


input and output in mass fractions 


* 1 input and output in mole fractions 


ICHUN 


Chemistry data units and unit conversion option. The GIM code chem- 
ically reacting flow model requires a variety of chemical information to 
describe the reactions and species involved. More often than not, these 
data as taken from JANAF tables, experimental calculations, and other 
sources is not in the correct set of units for use in the integrator module. 
Converting the data by hand before using it as input to the integration 
module can be tedious and time consuming. By selecting the appropriate 
ICHUN options and in conjunction with the conversion factors specified on 
card type C3, the user can input chemistry data in whatever units are con- 
venient and allow the program to perform all data conversion calculations. 


The initial value of ICHUN is 0 indicating no internal chemistry data 
conversion. The various ICHUN options can then be selected in any combina 
tion by adding the following values to ICHUN: 


ICHUN »0 no chemistry data conversion; chemistry data used 
as input to the program 

+ 1 general chemistry data conversion; appropriate 

data items are multiplied by user specified con- 
version factors (see card type C3) 

+ 2 per mole to per unit mass conversion; appropriate 
data items are divided by species molecular 
weight (Note: All GIM code chemically reacting 

flow model calculations are performed internally 
on a per unit mass basis. Therefore, all chem- 
istry data must be either input on that basis or 
converted to that basis. ) 


+ 4 reaction rate per particle to per mole conver- 
sion; appropriate reaction rate data are multi- 
plied by Avogadro's number 
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In addition to the value input for ICHUN, the sign of ICHUN also has 
significance. If ICHUN is input as a positive number, the chemistry data is 
printed out as input before any data conversion. If ICHUN is input as a 
negative number, the chemistry data are printed out both as input before 
data conversion and then again after data conversion has occurred. 

CARD TYPE C2 - Chemistry Control Tolerances 

FORMAT (4E10.0) 

Note: Card type C2 may be left blank if desired in which case the indicated 

nominal values will be used for CST0L(l-4). 

CST0L(1 ) (Nominal value = 1.0E-4) 

Trace species tolerance. Species whose mass fraction falls below the 
specified value of CSTOL(l) at any node in the flow field are considered to 
be trace species at that node. Trace species are included in all chemistry 
calculations but are not included in algorithmic decisions in the chemistry 
subroutines. 

CST0L(2 ) (Nominal value = 0.10) 

Explicit-to-implicit chemistry integration tolerance. If IXOPT is 
input as 0, the chemistry routines can switch from explicit integration of 
the species continuity equations to implicit integration. CST0L(2) controls 
the point at which the switch is made. After each explicit chemistry inte- 
gration step, the following criteria is examined: 

> CST0L(2) for all c £ > CSTOL(l) 
i-1, NSPEC 



If the criteria are satisfied, then the switch is made from explicit chem- 
istry integration to implicit chemistry integration on the next time step. 
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If the criteria are not satisfied, then the program continues to use ex- 
plicit chemistry integration. If IXOPT is not input as 0, then no switch- 
ing of chemistry integration scheme is possible and the value of CST0L(2) is 
irrelevant. 


CST0L(3) (Nominal value « 1.0E-10) 

Implicit-to-explicit chemistry integration tolerance. If IXOPT is 
input at 0, the chemistry routines can switch from implicit integration of 
the species continuity equations to explicit integration. CST0L(3) controls 
the point at which the switch is made. Before each implicit chemistry inte- 
gration step (immediately following the explicit provisional step), the 
following criteria are examined: 


MAX 



< CST0L(3) 


for all c i > CSTOL(l) 
i-1, NRSPC 


If the criteria are satisfied, then the switch is made from implicit chem- 
istry integration to explicit chemistry integration on the next time step. 

If the criteria are not satisfied, then the program continues to use im- 
plicit chemistry integration. Care must be exercised in selecting a value 
for CSTOL(3), particularly for equilibrium and nonequilibrium flows. If the 
value of CST0L(3) is too large, the program may switch from implicit chem- 
istry integration to explicit chemistry integration while the species pro- 
duction terms are still quite active. This will result in numerical insta- 
bility (time step criteria violation) and the calculation will "blow up." 

In flows with particularly violent reactions, a very small value for 
CST0L(3) is recommended in order to avoid this difficulty, e.g., 

CST0L(3) - 1.0E-10 

If IXOPT is not input as 0, then no switching of chemistry integration 
scheme if possible and the value of CST0L(3) is irrelevant. 
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CSTOL(A) (Nominal value =1.0) 


Chemistry iteration cutoff tolerance. During chemistry integration, 
the chemistry Eqs. (2.13) or (2.16) are solved iteratively using the method 
of successive substitutions or the Newton-Raphson method, respectively. 
These equations are iterated until the difference between successive 
iterates satisfies the following criteria: 


MAX 


u (£+l) _ U) 

i i 




< CST0L(4) 


for all c ± > CSTOL(l) 
i=l,NSPEC 


Limited experience has demonstrated that one iteration is usually sufficient 
to ensure numerical stability. Additional iterations are costly and not 
useful unless an accurate time transient is sought or unless the equations 
are stiff. Therefore, if only the steady state solution is desired and if 
the equations are not particularly stiff, CST0L(4) should be set to a large 
value (~1.0) so that the implicit chemistry routine will make only one 
iteration per time step. If CST0L(4) is set to a lower value, more itera- 
tions will be made on each time step but under no circumstances will more 
than 10 iterations be made 

CST0L(5) (Nominal value = 0.10) 

Chemistry Jacobian Matrix Update Tolerance. During chemistry integra- 
tion when Eqs. (2.13) or (2.16) are being solved iteratively, the Jacobian 
matrix (9H/9U) can be calculated only once at the beginning of the iteration 
process or recalculated on every iteration. The former approach is computa- 
tionally much less expensive than the latter but is not acceptable if the 
equations are stiff (9H/9U changing rapidly as a function of U). CSTOL(5) 
establishes the criteria to select which approach to use. If ISTlFFjtO, the 
following criteria is examined: 
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MIN 


(Convection + Diffusion Terms)^ 
(Species Production Terms)^ 


< CST0L(5) 

i=l ,NRSPC 


If the criteria are satisfied, the equations are assumed to be stiff and the 
Jacobian matrix is recalculated on every iteration. If the criteria is not 
satisfied or if ISTIFF=0, the equations are assumed to be non-stiff and the 
Jacobian matrix is calculated only once at the beginning of the iteration 
process. 


CARD TYPE C3 - Thermodynamic Conversion Factors 

FORMAT (4E10.0.4A8) 

Notes : 1. Input card type C3 only if MOD (CHUN ,2) ■ 0. 

2. If card type C3 is omitted (see it 1) or left blank, CPCF, HCF, 
SCF, and RRCF will be set to 1.0 internally and LCONV will be 
blanked . 

3. CPCF, HCF, and SCF sould not be used to convert from a per mole 
to a per mass basis. Use the ICHUN option for this purpose. 

CPCF 


Specific heat conversion factor. CPCF should be set equal to the 
factor required to convert specific heat data from input units to the proper 
integrator module units (ft 2 sec -2 °R _1 if IUNITS-1 or cm 2 sec -2 °K -1 if IUNITS 
* 2). Each of the specific heat ordinates input on card type C5 or read 
from file JANAF will be multiplied by CPCF before use in the program. 

Example : Specific heat data are taken from thermodynamic tables 

in units of cal/mole/°K and IUNITS « 2. 

CPCF should be set equal to 4.184 x 10 ^ g cm^ sec“^/cal. 


2-35 


HCF 


Enthalpy conversion factor. HCF should be set equal to the factor 
required to convert enthalpy data from input units to the proper integrator 
module units (ft 2 sec" 2 if IUNITS = 1 or cm 2 sec" 2 if IUNITS = 2). Each of the 
enthalpy ordinates input on card type 05 or read from file JANAF will be 
multiplied by HCF before use in the program. 

Example : Enthalpy data are taken from thermodynamic tables in units 

of kcal/mole and IUNITS = 2. 

HCF should be set equal to 4.184 x g cm 2 sec~^/kcal. 

SCF 


Entropy conversion factor. SCF should be set equal to the factor 
required to convert entropy data from input units to the proper integrator 
module units (ft 2 sec" 2 °R~ 1 if IUNITS = 1 or cm 2 sec" 2 0 K -1 if IUNITS = 
2). Each of the entropy ordinates input on card type C5 or read from file 
JANAF will be multiplied by SCF before use in the program. 

Example: Entropy data are taken from thermodynamic tables in units 

of BTU/mole/°R and IUNITS - 1. 

SCF should be set equal to 2.5036 x 10^ lbm ft^ sec~^/BTU. 


RRCF 


Activation energy conversion factor. RRCF should be set equal to the 
factor required to convert activation energy data (card type C9) from input 

1 or 


2 -2 

units to proper integrator module units (ft sec if IUNITS 
cm 2 sec" 2 if IUNITS 


2). Each of the reaction rate coefficients C. 
input on card type C9 will be multiplied by RRCF before use in the program. 


Example: Activation energy data are taken from thermodynamic tables 

in units of cal/mole and IUNITS = 1. 

RRCF should be set equal to 9.9287 x 10 1 lbm ft^ sec^/cal. 
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} 


LCONV 

Conversion factor label. LCONV is a 4 word label which is used on the 
printout as a mnemonic device to describe: the type of thermodynamic data 
conversions represented by the CPCF, HCF, SCF, and RRCF conversion factors. 
Any alphanumeric information can be used or LCONV can be left blasnk. LCONV 
is strictly a mnemonic device and in no way affects the calculations. 

CARD TYPE C4 - Species Identification Label, Molecular Weight, and Heat 
of Formation 

FORMAT (A6,4X,2E20.0) 

Notes : 

1. Input card types C4 and C5 only if NTIP > 0 on card Cl. 

2. For NTIP > 0, the input sequence is one card type C4 followed by 
several cards of type C5; then repeat the sequence NSPEC times, 
card C4 - then cards C5, card C4 - then cards C5, etc. 

3. The number of cards of type C5 to be input after each card type 
C4 is equal to (NTIP-l)/2 + 1. 

4. In the subsequent input description, the subscript I runs from 1 
to NSPEC, indicating the sequence of NSPEC groups of cards type 
C4 and C5. 

LSPID(I) 

Species identification label for the I s species. Any alphanumeric 
data up to six characters in length may be used. The explicit chemical 
formula for the species may be used but this is not required. For example, 
methane might be represented by either CH4 or MTHANE. Note however that 
whatever label is used for the 1^ species in LSPID(I) must also be used 
on card type C8 to represent that species in any reaction mechanisms. 
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WM(I) 


Molecular weight of the I t ^ 1 species. 
HFM(I) 


Heat of formation of the I^ 1 species. HFM(I) should be input in 
units consistent with the enthalpy ordinates input on card type C5. Like the 
enthalpy ordinates, HFM(I) will be multiplied by the HCF conversion factor 
if M0D(ICHUN,2) ^ 0 and will be divided by WM(I) if M0D(ICHUN/2 , 2) t 0. 

CARD TYPE C5 - Species Thermodynamic Properties Data 

FORMAT (8E10.0) 


Notes 

1. The number of cards of type C5 to be input (after each of the 
NSPEC card type C4) is equal to (NTIP-l)/2 + 1. 

2. In the subsequent input description, the subscript J runs from 1 
to NTIP (incremented by 2) , indicating the sequence of NTIP 
thermodynamic abscissa and ordinate points for each of the NSPEC 
species. 


TIP(J) 


Temperature abscissa (J th value). Note that TIP is not subscripted 
with the subscript I which implies that the same set of TIP must be input 
for each of the NSPEC species on card type C5. The NTIP values of TIP 
should span the range of temperatures expected over the entire flow field 
during the time integration. If the temperature at any node in the flow 
field exceeds the range of the TIP values during integration, an execution 
diagnostic message is printed and the calculation terminates (no extrapola- 
tion is performed). TIP values need not begin at zero degrees nor be 
equally spaced on input. However, TIP and its associated data ordinates 
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will be extrapolated to zero degrees and interpolated to equally spaced 
abscissas internally before time integration begins. TIP should be input in 
unit of °R (IUN1TS = 1) or °K (IUNITS - 2). 

CPHS( J > K > I) t K~l,3 


Species thermodynamic data ordinates corresponding to the J ^ tem- 
perature abscissa, TIP(J), for the I c species. CPHS(J,K,I) data consist 
of the following items: 


CPHS(J ,1,1) «* specific heat at constant pressure 


W 


CPHS(J,2,I) = sensible enthalpy * _ / ^ c (£) d£ 

r J p . 

o *i 

CPHS(J,3,I) = entropy - s i< T j) 


The CPHS thermodynamic data ordinates for all species should be input in 
consistent units and are subject to data conversion via the CPCF, HCF, and 
SCF factors if M0D(ICHUN,2) / 0 as well as per mole to per unit mass conver- 
sion if MOD (I CHUN/ 2, 2) + 0. 

TIP( J+l 


Temperature abscissa (J+l 


th 


value) . 


See TIP(J) above. 


CPHS ( J+l, K, I) , K-1,3 


Species thermodynamic data ordinates (J+l 
above . 


th 


values) . 


See CPHS(J.K.I) 
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CARD TYPE C6 - Catalytic Species Identification Label 

FORMAT (A6) 


Notes 

1. Input card types C6 and C7 only if NCATS >0 and NREAC > 0 on 
card Cl. 

2. For NCATS > 0, the input sequence is one card type C6 followed by 
jne or more cards of type C7; then repeat the sequence NCATS 
times, card C6 - then cards C7, - then cards C7, etc. 

3. The number of cards of type C7 to be input after each card type 
C6 is equal to (NCATS-1)/16 + 1. 

4. In the subsequent input description, the subscript I runs from 1 
to NCATS, indicating the sequence of NCATS groups of cards type 
C6 and C7. 


LSPID(I-fNSPEC) 

Species identification label for the I c catalytic species. Any 
alphanumeric data up to six characters in length may be used. The explicit 
chemical formula (if any) may be used, but this is not required. Note that 
whatever label is used for the I C ^ catalytic species in LSPID(I+NSPEC) 
must also be used on card type C8 to represent that catalytic species in any 
direction mechanism. 

CARD TYPE C7 - Catalytic Species Weighting Factors 

FORMAT (16F5.0) 

Notes 

1. The number of cards of type C7 to be input (after each of the 

NCATS card type C6) is equal to (NCATS-1)/16 + 1. 

2. In the subsequent input description, the subscript J runs from 1 

to NSPEC, indicating the NSPEC weighting factors for each of the 

NCATS catalytic species. 
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WFM(I.J) 


Weighting factor for reactive species J associated with the cata- 
lytic species. During reaction rate calculations, the presence of catalytic 
or ’’third-body" species is included as a fictitious mass fraction consisting 
of the weighted sum of the reactive species mass fractions. At each node in 
the flow field the 1^ catalytic species mass fraction is determined as 
follows: 

NSPEC 

C cat. ^ M ij * c reac, 

1 , j 


No restrictions apply to the values of WFM(I,J). 

CARD TYPE C8 - Reaction Mechanism Information 

FORMAT (211, 78A1) 


Notes 

1. Input card type C8 only if NREAC > 0 on card Cl. 

2. For NREAC > 0, the input sequence is one card type C8 for each of 
the NREAC reactions under consideration. 

3. In the subsequent input description, the subscript J runs from 1 
to NREAC, indicating the sequence of NREAC cards of type C8. 

RTYPE(I.J), 1-1,2) (Type INTEGER variable) 

Reaction rate constant type associated with reaction J. RTYPE(1,J) is 
associated with the forward reaction rate constant while RTYPE(2,J) is asso- 
ciated with the backward reaction rate constat. RTYPE values are selected 
depending on the algebraic form of the corresponding reaction rate constant 
calculation, as indicated below: 


2-41 


RTYPE = 0 
= 1 
- 2 

- 3 

- 4 

- 5 

- 6 


k = 0.0 (no reaction - backward reaction only) 

k = A (constant) 

k = A(T-T q ) B 

k = A exp(C/!R/(T-T 0 )) 

k = A(T-T q ) B exp(C/K/(T-T o )) 

k = A(T-T q ) B exp(C/5l/(T-T o ) D ) 

k « calculated from equilibrium constant (backward 
reaction only (Eqs. (2.4e) through (2.4g)), 


Note that the value of RTYPE ( I, J) must be in the following range: 

1 < RTYPE(1,J) < 5 
0 £ RTYPE ( 2 , J ) £ 6 

LREAC(I), 1-1,78 


Reaction mechanism equation for reaction J. Each of the NREAC reac- 
tions should be input via LREAC on separate cards of type C8. Each equation 
can consist of up to 78 characters including blanks and must adhere to the 
following format: 

SUM OF REACTANTS « SUM OF PRODUCTS 

where each reactant and each product consists of an optional integer coef- 
ficient and a LSPID species identification label. If the optional integer 
coefficient is present, it must be separated from the LSPID labvel by an 
asterisk (*). Individual reactants and products are separated by a plus 
sign (+)• In this format each reaction equation should resemble a FORTRAN- 
like expression, e.g., 


H2 + 02 = OH + OH 
2 * H20 « 2 * OH + H2 


As explained above, the characters *, +, and « have special significance and 
should only be used to delimit coefficient and species (*), different 
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species (+) , and reactants and products (=). In particular, these three 
special characters should not be used in any fashion in either the reactive 
species identification label (LSPID) on card type C4 or the catalytic 
species identification label on card type C6. Blanks are not significant 
and may be used for spacing as desired. 

CARD TYPE C9 - Reaction Rate Coefficients 

FORMAT (5E10.0) 

1. Input card type C9 only if NREAC > 0 on card type Cl. 

2. For NREAC > 0, the input sequence consists of one or two cards of 
type C9 for each of the NREAC reactions under consideration. One 
card is always input; the second card is input only if RTYPE(2,J) 
is greater than 0 and less than 6 on card type C8, 

3. In the subsequent input description, the subscript J runs from 1 
to NREAC. 

RRATC(I,1,J) 1=1,5) 


Forward reaction rate coefficients. RRATC(I,1,J) are the reaction rate 
coefficients for the forward reaction of reaction mechanism J. The 
RRATC(I , 1 , J) coefficients correspond to the rate coefficients in Eq. (2.4c) 
with 


RRATC(1,1, J) = Aj 
RRATC(2,1, J) = Bj 
RRATC(3,1,J) = CJ 
RRATC(4,1,J) - Dj 
RRATC(5,1, J) = r 


pre-exponential factor 
temperature exponent 
activation energy 
temperature exponent 
reference temperature 


All five of these coefficients are input for each of the NREAC reactions. 
The value of RTYPE(1,J) on card type C8 determines which of the five are 
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significant in the calculation of the forward reaction rate. RRATC(1,1,J) 
is subject to per particle to per mole conversion if M0D(ICHUN/4, 2) = 0 and 
RRATC(3,1,J) is subject to data conversion via the RRCF factor if 
M0D(ICHUN,2) = 0. 

RRATC(I j 1_ l J) 1=1,5 


Backward reaction rate coefficients. RRATC(I,2,J) are the reaction 
rate coefficients for the backward reaction of reaction mechanism J. The 
RRATC(I,2,J) coefficients correspond to the rate coefficients in Eq. (2.4d) 
with 


RRATC(1 , 2 , J) = Aj' 
RRATC(2 , 2 , J) = B'j 
RRATC(3,2,J) = Cj 
RRATC(4 ,2, J) = Dj' 
RRATC(5 , 2 , J) = 


pre-exponential factor 
temperature exponent 
activation energy 
temperature exponent 
reference temperature 


All five of these coefficients are input for each of the NREAC reactions 
when RTYPE(2,J) is greater than 0 and less than 6. Otherwise, this second 
card type C9 is omitted and the backward reaction rate is either identically 
zero (RTYPE(2,J) - 0) or is calculated from the equilibrium constant 
( RTYPE ( 2 , J ) - 6) using Eqs. (2.4e) through (2.4g). If RTYPE(2,J) is greater 
than 0 and less than 6, the value of RTYPE(2,J) on card type C8 determines 
which of the RRATC(I,2,J) coefficients are significant in the calculation of 
the backward reaction rate. RRATC(1,2,J) is subject to per particle to per 
mole conversion if M0D(ICHUN/4, 2) 4 0 and RRATC(3,1,J) is subject to data 
conversion via the RRCF factor if M0D(ICHUN,2) 4 0. 
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CARD TYPE CIO - Species Initialization Control 

FORMAT (415) 


Notes 

1. Any number of cards type CIO may be used to initialize the 
species mass/mole fraction distribution. All nodes can be input 
on a single card or each node can be input on a separate card. 

The usual case is somewhere between these two extremes. 

2. If I TYPE * 0 is input on card type CIO, then (NSPEC-l)/8 + 1 
cards of type Cll must immediately follow. If ITYPE f 0, then no 
cards of type Cll are input. 

3. A -1 card (columns 4-5) must be input as the last card in a type 
CIO sequence to terminate reading of species mass/mole fraction 
initial data. 

4. Cards type CIO (and Cll if required) are not input for a restart 
case (ISTART f 0) unless changes are being made. The -1 card 
must be presnt even on a restart case. 

5. For a quasi-parabolic run (METHOD > 2), one sequence of cards 
type CIO (and cll if required) must be input with the other chem- 
istry data (cards Cl through C9) to initialize the first plane. 

If there are added zones in the problem, then additional se- 
quences of cards type CIO (and Cll if required) may be needed to 
Initialize species mass/mole fraction in the added zone plane. 

If this is the case, additional cards of type CIO (and Cll if re- 
quired) should follow after Integrator module cards type 16 and 
16a. 

6. If IMOLE - 0 on card type Cl, then species mass fractions must be 
input on card type Cll. If IMOLE - 1 on card type Cl, then 
species mole fractions must be input on card type Cll. 


NJ 


Node number of the first nodal point to be initialized by this card 


type CIO. (NJ ■ -1 terminates the input of card CIO). 
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INC 


Node number increment to NJ to be used for inputting a sequence of 
nodes on one card. Set INC = 0 if only one node is to be set by this card. 
If INC f 0, it will be added to NJ to determine the next node to be initial- 
ized . 


NTOT 


The total number of nodes to be set by this card type CIO. 


ITYPE 


Indicates the type of species mass/mole fraction initialization is to 
be done. 

ITYPE = 0 allows the user to input the initial species mass/ 
mole fraction distribution on cards of type Cll 

* 1 sets the species mass/mole fractions for this card 

type CIO input to the same values as input on the 
last card type Cll 

* 2 allows the user to code and execute subroutine USERIP 

in order to initialize the species mass/mole fractions. 


CARD TYPE Cll - Species Mass/Mole Fraction Initial Values 

FORMAT (8E10.0) 


1. Cards of type Cll are input following a card type CIO only if 
ITYPE ■ 0 is input on card type CIO. 

2. If ITYPE * 0 on card type CIO, then (NSPEC-l)/8 + 1 cards of type 
Cll must follow. 

3. Cards of type Cll are not input on a restart case unless changes 
are being made to the species mass/mole fraction distributions. 
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1 

4. If IMOLE ■ 0 on card type Cl, then species mass fractions must be 
input on card type Cll. If IMOLE * 1 on card type Cl, then 

' species mole fractions must be input on card type Cll. 

5. In the subsequent input description, the subscript I runs from 1 
to NSPEC, indicating the NSPEC species mass/mole fraction initial 
values associated with the nodes specified on card type CIO. 


CS(I) 


Initial mass/mole fraction value for species I at those noses specified 
on preceedlng card type CIO. CS(I) must be in the following range: 

0.0 _< CS(I) _< 1.0 

In addition, the CS(I) must obey the following relationship: 


NSPEC 

^ ^C 1 ) " 1 *° 
1-1 


2.6 SAMPLE CALCULATIONS 


The results from two sample calculations using the GIM code chemically 
reacting flow model are presented in this section. These two test cases 
include: 


1. Premixed hydrogen-air ignition in a one-dimensional duct-like 
geometry using a five-species, twoi-reactlon global chemistry 
model developed for this problem by the method of Appendix A; 
and 

2. The hydrogen-air parallel-injection case of Burrows and Kirkov 
(Ref. 2-5) using the five-species, two-reaction global chem- 
istry model developed by Rogers and Cbinitz (Ref. 2-6). 


The first case consists of the ignition of premixed hydrogen and air in 
a one-dimensional 40 cm x 10 cm duct-like configuration. Figure 2-3 shows 


I 


2-47 



Fig. 2-3 - Premixed I^-Air Case (Computational Grid) 


the geometry and computational mesh used in this calculation. The global 
reaction mechanism was synthesized using the method of Appendix A for this 
geometry and flow conditions. The reaction mechanism consists of the five 
species N2, 02, H2, H20, and OH and the following two reactions 


H2 + 02 ZZT OH + OH 


k_ = 4.171 x 10 7 (T - 1000. 0) 0,98 
F i 


H20 + H20 — ~0H + OH + H2 


k Fi = 1.054 x 10 2 T 5 * 87 exp[- 


3 

where kp and kp are in units of cm /mole/sec. The inflow boundary 
conditions are given below: 

P = 3.204 x 10 4 gm/cm 8 

u ■ 1.464 x 10 8 cm/sec 

T - 1000 K 

6 2 

P * 1.013 x 10° dyne/cni * 1 atm 

M = 2.25 
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0.76149 


23123 


0.00728 


0 


C 


C 


02 


C 


H2 


C H2 " C 0H 


These same conditions were also used as the initialization for the remainder 
of the flow field. The integration was performed using the explicit chem- 
istry solution option until the OH concentration reached a sufficient level 
for ignition to begin. At that point the code automatically switched to the 
implicit chemistry solution option and the calculation continued. The time 
step for this calculation was maintained at 0,95 times the flowfield CFL. 

The results from this calculation presented in Figs. 2-4 through 2-8 
show a sharp ignition front where the hydrogen-air combustion occurs. Fig- 
ure 2-8 also presents a comparison with the results obtained from a spatial 
marching chemistry code (ALFA) for the same calculation. The discrepancy 
between the GIM code results and the results from the ALFA code solution are 
attributable to the inadvertent use of slightly different reaction rate 
coefficients in the two calculations. 

The second sample calculation consists of the parallel injection case 
of Burrows and Kurkov (Ref. 2-5). In this case a cold sonic stream of 
hydrogen is injected parallel to a supersonic freestream of hot air. The 
resulting mixing and combustion are calculated and compared with the experi- 
mental data. Figure 2-9 shows the geometry of the test section used by 
Burrows and Kurkov to obtain the experimental results. A sparse version of 
the computational grid used in this calculation is shown in Fig. 2-10. This 
grid only displays every third node in the x-direction and every fifth node 
in the y-direction. The actual gjrid used in the calculation is thus much 
more dense, consisting of 109 cross planes with 56 nodes per plane. The 
nodes are clustered around the hydrogen jet near the inflow region in order 
to resolve the mixing that takes place. A Baldwin-Lomax turbulence model 
(see Section 5) was used to simulate the mixing process. Because the mixing 
and combustion takes place in a region close to the lower channel wall, the 
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X (cm) 


Fig- 2-7 Premixed I^-Air Case (Pressure vs X) 
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Fig. 2-8 - Premixed l^-Air Case (Log Mass Fraction vs X) 


i 
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I 


Test section intermediate 
measurement station, 

Test section initial 18. 3 cm Test section exit 



Fig. 2-9 - Parallel Injection Case (Test Section Geometry) 



Fig. 2-10 - Parallel Injection Case (Computational Grid (Coarse)) 
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upper boundary was treated as an inviscid f ree-slip/tangency boundary in 
order to avoid the necessity of resolving the boundary layer there. The 
lower channel wall was treated at a no-slip constant temperature boundary. 

The flowfield conditions for this case are given below: 


Hydrogen Jet 


Free Stream 


Mach Number 
Temperature (K) 
Velocity (cm/sec) 

2 

Pressure (dyne/cm ) 
Mass Fractions: 


1.00 

254 

1.216 x 10 5 
1.1 x 10 6 


C 

C 

C 

C 

C 


H2 

02 

N2 

H20 

OH 


1.000 

0 

0 

0 

0 


2.44 

1270 

1.764 x 10 5 
1.1 x 10 6 


0 

0.258 

0.486 

0.256 

0 


The global hydrogen-air reaction mechanism developed by Rogers and 
Chinitz (Ref. 2-6) was used in this second sample calculation. This partic- 
ular reaction mechanism consists of the five species N2, 02, H2, H20, and OH 
and the following two reactions: 

H2 + 02 ^ 2 OH k = 1.3823 x 10 48 T _1 ° exp[ ~ 4885 ] 

RT 

k_ 

F 

2 0H + H2 — l 2 H20 k_ = 2.7166 x 10 64 T~ 13 exp[ ~ - 4 j 3Q0 ] 

* r 0 RT 
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3 6 2 

where k is in units of cm /mole/sec and k is in units of cm /moleVsec. 

F 1 F 2 
These reaction rate coefficients resulted in a set of very stiff species 

continuity equations and required the use of the implicit chemistry solution 

option (1X0PT = -1) and the stiff equation option (ISTIFF = 1) throughout 

the calculation. The integration was performed using the flowfield results 

from a standard GIM code binary gas integration run (with the same geometry 

and boundary conditions) as te initialization for this reacting case. The 

flow field was also initially "seeded” with small concentrations of OH in 

order to lessen the stiffness of the species continuity equations. The time 

step for this calculation was maintained at 0.75 times the flowfield CFL. 


The results from this second calculation are presented in Figs. 2-11 
through 2-16 for the test section exit station, x 55 35.6 cm. The results 
show a very hot region of hydrogen-air combustion approximately 2.0 cm above 
the lower channel wall. Figure 2-16 shows a comparison of the GIM code 
solution with the Burrows and Kirkov experimental data. The comparison is 
relatively good although some of the calculated peaks in the H20, N2, and 02 
distributions are not quite as high as those indicated by the experimental 
data. 
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Fig. 2-11 - Parallel Injection Case (Density vs Y) 
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Fig. 2-12 - Parallel Injection Case (Velocity vs Y) 
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ITER NO. 7000 



Fig. 2-13 - Parallel Injection Case (Temperature vs Y) 
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Fig. 2-14 - Parallel Injection Case (Pressure vs Y) 
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Fig. 2-15 - Parallel Injection Case (Log Mass Fraction vs Y) 
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2 . 7 NOMENCLATURE 


English Symbols 


a acoustic velocity 

A temperature/enthalpy quadratic equation coeffi- 

cient; Eq. (2.21c) 

Aj forward reaction rate coefficient, reaction j; Eq. 

(2.4c) 

A 1 :' backward reaction rate coefficient, reaction j; Eq. 

(2.4d) 

B temperature/enthalpy quadratic equation coeffi- 

cient; Eq. (2.21d) 


B 


j 


forward reaction rate exponent, reaction j; Eq. 
(2.4c) 


B 


»» 

j 


backward reaction rate exponent, reaction j; Eq. 
(2.4d) 


c i 


mass fraction of species i 


Cp specific heat at constant pressure of species i 

C temperature/enthalpy quadratic equation coeffi- 

cient; Eq. (2.21e) 


forward reaction rate coefficient, reaction j; Eq. 
(2.4c) 


C 


ii 

j 


backward reaction rate coefficient, reaction j; Eq. 
(2.4d) 


C: 


denotes corrector step 


CH: 


denotes chemistry step 




forward reaction rate exponent, reaction j; Eq. 
(2.4c) 


Dj backward reaction rate exponent, reaction j; Eq. 

(2.4d) 

artificial numerical diffusion coefficients for 
1 2 3 xj_, X2, and X3 components 
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JS £ 


exp 

e j 

€ 

AG i 

h f. 

1 

hi(T) 


real diffusion coefficient 

sum of real and artificial diffusion coefficients 

exponential function 

vector of flux terms; Eq. (2.6b) 

total internal energy per unit mass 

change in Gibbs' free energy for reaction j; Eq. 
(2.4g) 

enthalpy of formation per unit mass of species i 
enthalpy per unit mass of species i at temperature T 


h i,j 



enthalpy per unit mass ordinate of species i at 
temperature abscissa Tj 

enthalpy per unit mass at reference temperature 
T 0 of species i 


hf enthalpy of formation per mole of species i 

h Q enthalpy per mole at reference temperature T 0 of 

species i 


(hi - h 0i ) 


enthalpy per mole based on reference enthalpy h 0 
of species i 


kflow field 


enthalpy per unit mass as computed from flowfield 
variables; Eq. (2.18a) 


h 


specie 


enthalpy per unit mass as computed from species 
variables; Eq. (2.18b) 


H vector of production terms, Eq. (2.6b) 

(3H/3U) Jacobian matrix of vector H 


I 


identity matrix 


J 


number of discrete equally spaced temperature 
abscissa 


k 


coefficient of thermal conductivity 



backward reaction rate constant for reaction j; Eq. 
(2.4d) 
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♦r-> 

<4-1 

* 

forward reaction rate constant for reaction j; Eq. 
(2.4c) 


concentration equilibrium constant for reaction j; 
Eq. (2.4e) 

s 

pressure equilibrium constant for reaction j; Eq. 
(2.4f) 

m i 

molecular weight of species i 

M 

number of reactions 

M 

implicit coefficient matrix; Eq. (2.23b) 

N 

number of species 

P 

pressure 

P: 

denotes predictor step 

q 

velocity magnitude 


heat flux per unit area for x^, X£, X3 co- 
ordinate directions 

R 

N ^ 

gas constant for mixture ■ J c — 

i.-l 1 m i 

51 

universal gas constant 

RHS 

implicit right-hand side vector; Eq. (2.23b) 

s"i 

entropy per mole of species i 

t 

time 

At 

time step 

T 

temperature 

T j 

temperature abscissa, j 1 * 1 value 

T ° 

reference temperature 

t' 

J 

forward reaction rate reference temperature, re- 
action j; Eq. (2.4c) 

—II 

°j 

backward reaction rate reference temperature, re- 
action j, Eq. ( 2 . 4d ) 

AT 

temperature increment 

! 
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u l, u 2» u 3 

U 

'U 

6u n+1 


AU nH ' 1 


Wi 

x l> x 2> x 3 

X j 

Greek Symbols 



e 


X 


P 

V 


V 


ij 


A v 


ij 


T 


ij 


Operators 


B[ ] 

3 B /3xj 


3 F /9xj 


velocity components 

vector of conserved variables; Eq. (2.6b) 

implicit time derivative vector; Eq. (2.23b) 

finite difference time derivative of conserved 
variables at time level n+1 

change in conserved variables at time level n+1 
chemistry production term for species i; Eq. (2.4b) 
coordinate directions 

reaction components defied by Eq. (2.25b) 

Kronecker delta = 1 for i = j, =0 for i *= j 
artificial diffusion constant 
coefficient of bulk viscosity 
mass density 
coefficient of viscosity 

stoichiometric coefficient for species i as a re- 
actant in reaction j 

stoichiometric coefficient for species i as a pro- 
duct in reaction j 

it i 

Vij - Vij 

stress tensor; Eq. (2.2b) 
boundary condition operator 

backward spatial difference operator for coordinate 
direction j 

forward spatial difference operator for coordinate 
direction j 
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Superscripts 


(£) 

implicit chemistry 
n+1 time 

iteration for time level 

n 

level n 


n+I 

first provisional time 

level (predictor step) 

n+1 

second provisional time 

level (corrector step) 

n+1 

time level n+1 


overbar 

molar basis 
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3. THE GIM HYPERBOLIC STEADY STATE EULER SOLVER 


3 . 1 METHODOLOGY 

The success of the General Interpolants methodology in developing com- 
puter codes for three-dimensional flowfield calculations in arbitrary geo- 
metric domains has lead to an elliptic unsteady Navier-Stokes solver for 
complete viscous, compressible flow analysis and an iterative quasi-unsteady 
parabolic Navier-Stokes solver for viscous flows with a predominant flow 
direction along which second order effects can be neglected. Although the 
latter method has been used extensively (Refs. 3-1 through 3-4) to solve the 
completely inviscid supersonic Euler equations, the iterative procedure is 
unnecessary for such flows since the equations are hyperbolic in the flow 
direction and it adds considerable computational cost. The impetus for using 
the quasi-parabolic code in the inviscid mode was to take advantage oi the 
geometric versatility of the GIM approach. This section reports on a 
non-iterative hyperbolic steady-state Euler equation solver developed through 
the General Interpolants Method which has the following among its many 
feature s : 

• The geometric versatility of all GIM derivatives 

• The second order accuracy, efficiency and highly vectorizable char- 
acteristics of a MacCormack predictor-corrector solution algorithm 

• The advantages of the new technique for progressive assembly of the 
GIM difference analogs (see Appendix D) 

• Continuous monitoring of the marching stepsize with the capability 
to interpolate intermediate cross planes to prevent violation of 
stability limits 
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• Entropy consistent inviscid boundary conditions 

• Shock-capturing techniques 

• Integration of the equations in the Cartesian physical domain. 

This solver has proven to be quite accurate and efficient in axisym- 
metric, two- and three-dimensional test cases. The next subsection details 
the finite difference model. The following two subsections describe the 
stepsize control mechanism and the boundary conditions, respectively. 

3.1.1 The Finite Difference Model 


The Euler equations for three-dimensional compressible inviscid gas 
dynamics are 


where 


3U . 3E , 3f , 3G 
"ST + "5x + ITy + Tz 
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and where 


p = mass density 

H = total enthalpy 

P *■ pressure 

£ = total energy 

u = x-component of velocity 

w = z-component of velocity 

v * y-component of velocity 

x,y,z,t = space and time coordinates 


and 


c 1 -^) p 


H - 


u 2 + v 2 + w 2 


(ideal gas law) 


Application of the General Interpolants methodology (Ref. 3-5) using 
the weight functions derived for the progressive assembly of the derivative- 
taking analogs (Appendix D) results in the following finite difference 
analog of the Euler equations at node N. 


A(x,y,z) 
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where the terms in brackets are the assembled finite difference Jacobians: 


and 
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is the element finite difference Jacobian for element e as described in 
Appendix D. The are the arbitrary weighting factors which determine the 
direction of the differencing and 



(3.3) 


One choice of weighting factors which satisfy Eq. 
tri-linear interpolants: 


(3.3) is a set of 
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where the 0^ govern the differencing direction in as follows: 

0^ = 0 - forward differences 

0^ = 1/2 - centered r|^ difference? 

0^ = 1 - backward differences 

Thus, the coefficient of Ujj in Eq. (3.2) can be viewed as an interpolated 
control volume for the discrete difference model. 


For the steady-state Euler equations we set = 0. The finite dif- 
ference model is a spatial marching version of the MacCormack explicit 
predictor-corrector scheme (Ref. 3-6). To implement this scheme we choose 
0^ = 0 for both steps and 02 and 9g are chosen to yield the classic 
alternating differences on the predictor and corrector steps. The choice of 
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cross-flow differencing is completely arbitrary and user input, 
ing algorithm may be written 


The result- 


Predictor : 
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(3.5a) 


Corrector: 
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(3.5b) 


where, again, the terms in brackets are the two and three-dimensional 
Jacobian determinants and 6^ is the appropriate two-point difference oper- 
ator in the direction for element e. a & and a & are the predictor 
and corrector weighting factors, respectively. The overbar (— ) represents 
predicted values. The marching proceeds down ri^ grid lines from station n 
to station n+1. 


Notice that the metric information is evaluated at station n on both 
steps. This is done, first, because the advanced geometry lies outside the 
domain of dependence of both the partial differential equations and the 
difference equations, and, secondly, from the realization that the actual 
advancement of the solution occurs on the corrector step. The primitive 
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variables are decoded from the flux vector E assuming supersonic flow in the 
x direction* 

The geometric versatility of the other GIM derivatives was maintained 
by retaining the terms in 6^F and S^G. This does not require the cross- 
flow mesh to be orthogonal to the x axis* It does, however, require an 
iterative procedure for non-orthogonal meshes. This iteration converges 
quickly ( <5 iterations) and is a small penalty to pay for the increased 
geometric power* 


This algorithm readily lends itself to vectorization on the CYBER 200 
series of supercomputer. This efficiency along with the proven accuracy of 
the technique far outweigh the limitations placed on the marching stepsize 
by its explicit nature and preclude the use of much less computationally 
efficient implicit methods. 

3.1,2 Hyperbolic Stability Analysis 

As is the case in all explicit hyperbolic marching algorithms, the mag- 
nitude of the spatial or marching step size used in the GIM code steady 
state marching algorithm is constrained by the hyperbolic nature of the flow 
field. The maximum allowable spatial step size is given by the GFL cri- 
teria, i.e., the domain of dependence of the partial differential equations, 
which comprises the area enclosed by the intersection of the Mach conoid 
through the solution point and the previous initial value surface, must be 
contained within the domain of dependence of the finite difference equations 
(Ref. 3-7). This concept is shown geometrically in Fig. 3-1. 


The equation of the Mach conoid is given by (Ref. 3-7); 


222 2222 2229 9 

[u - (V - a )] (dx) + [v - (V - a )] (dy) + fw - (V - a )] (dz) 


+ 2 uv dx dy + 2 uw dx dz + 2 vw dy dz « 0 


(3.6) 
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where V is the velocity magnitude and a is the speed of sound. Let the 
intersection of the Mach conoid from a node on the solution surface (point 1 
in Fig, 3-1) with the grid line passing through the previous solution node 
and any connected node (point 0 and points 3,4,5, and 6 in Fig, 3-1) be 
designated by x and let its position be defined in terms of a parameter 9 , 
then 
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where subscript c denotes any connected node. On the Mach conoid, let the 
differentials be approximated by differences, i.e., 
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Then Eqs. (3,7) and (3.8) together yield 
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Equation (3.6) can then be approximated by 


AAx^ + BAxAy + C£y +DAyAz + EAz +FAxAz 


(3.10) 


where 
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If the definitions of4x,dy, and 4z given by Eq. (3.9) are substituted into 
Eq. (3.10) and like powers of 6 are collected, the results is given by 


a e 2 - BQ + d - 0 

where 

a = A Ax 2 + B Ax Ay + C Ay 2 + DAy. Az + EAz 2 + FAx Az 
c c c c cc c cc 

y5 • 2AAxjAx fi + B(Ax^ Ay c + Ax c Ayp + 2C Ay^ Ay fi 

+ D(AyjAz c + Ay c Azj) + 2E AzjAz £ + F(AXjAz c + Ax £ Az^) 

<3- ■ AAx 2 + B AXjAy^ + CAy^ 2 + DAyjAz^ + EAz^ 2 + FAx^ Az^ 


(3.11) 


Equation (3.11) is a simple quadratic equation in 9 with the solution given 
by 
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(3.12) 
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If 0 <_ 1, then the differential domain of dependence is contained within 

the difference domain of dependence (Fig. 3-la) and the step between node 0 
and node 1 is stable. If 0 >1, then the differential domain of 

dependence is not entirely contained within the difference domain of 
dependence (Fig. 3-lb) and the step between node 0 and node 1 will be 
unstable. Because the above analysis is approximate, a user supplied safety 
factor is applied to the value of 9 determined by Eq. (3.12) in order to 
ensure that marginal stability does not present a problem. 

If the above analysis indicates that the marching step size associated 
with the previously computed "primary** geometry is unstable, the GIM code 
steady state marching algorithm can automatically supply intermediate or 
"secondary" solution surfaces to ensure that the spatial marching proceeds in 

a stable fashion. The geometry of any required secondary solution surfaces 
is obtained via linear interpolation between primary geometry surfaces. 

3.1.3 Decode Procedures 


After the E vector has been obtained at a new station, the primitive 
variables p, u, v, w, H, and P must be decoded. This is accomplished in the 
following six steps: 
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(3.13) 
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then the primitive variables are given by 


1. v 

2. w 

3. H 

4. u 

5. P 

6. P 


V E i 

V E i 

V E 1 


E /u 


P (-^-) [H - y (u‘ + v‘ + w-)] 


Y ^2 

y+t ^ it ^ 


+ V ( rlr )2 4r> - 


2(^-) [H - j (v 2 + w 2 )] (3.14) 


This decode procedure is used after both the predictor and corrector steps 
in order to obtain the primitive variables from the E vector. 

3.1.4 Boundary Conditions 


After the primitive variables are obtained from the decode procedure, 
the boundary conditions must be enforced. In the inviscid flow of an ideal 
gas this consists of requiring that the component of velocity normal to any 
solid boundary vanish, l.e., the so-called free-sllp tangency boundary con- 
dition. In the GIM code steady state marching algorithm, this is accom- 
plished using a variation of Abbett's method (Ref. 3-8) described below. 

At any solid boundary, the normal component of velocity must vanish, 

l.e. , 


q • n ■ 0 


(3.15) 
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where q is the velocity vector with Cartesian components (u,v,w) and n is 
the outward pointing (into the computational domain) unit normal vector with 
components (n x , n^, n z ). In general, the velocity vector components 
obtained from the decode procedure, Eq. (3.14), will not satisfy Eq. (3.15) 
identically but will be rotated out of the surface tangent plane by a small 
angle S. This angle can be determined by 

6 38 sin * (q • n /q) (3.16) 


where q is the magnitude of the velocity. 

If <5 is positive, then Abbett's method indicates that an expansion wave 
is necessary to rotate the velocity vector back into the tangent plane but 
if 6 is negative, then a compression wave is required. In practice, this 
means correcting the values of p, u, v, w, and Pj calculated from 
MacCormack’s method, by using shock wave and expansion wave properties. 

This is accomplished as follows: 


1. Calculate p, u, v, w, H, and P from Eqs. (3.5a, 3.5b) and (3.14). 

2. Calculate 6 using Eq. (3.16). 

3. If 6 is positive, go to step 5. 

4. $ negative indicates compression wave required. 

a. Compute oblique shock wave angle 0 for deflection angle |fi | and 
Mach number ■ q^/(YP/p) using Eqs. (150) from Ref. 3-9: 
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There are three roots of Eq. (3,17), the smallest of which cor- 
responds to a decrease in entropy and is therefore discarded. 
The largest root corresponds to a strong oblique shock and is 
also discarded. The intermediate root corresponding to a weak 
oblique shock is the one of interest. 

b. Compute the corrected static pressure P c using Eq. (128) from 
Ref. 3-9: 


p = p 2 YM 2 sin 2 9 - (Y~l) 
c Y+l 


(3.18) 


c. Compute the entropy Increase across the shock using Eq, (144) 
from Ref, 3-9 : 



2YM 2 sin 2 9 - (Y-l) 

-Y In 

(Y+1)M 2 sin 2 0 

L Y+l 

_(Y-1)M 2 sin 2 0 + 2 _ 



(3.19) 


d. Compute the corrected density p c using fundamental 
relationships: 


Pc 



1/Y 



(3.20) 


where subscript u indicates values immediately upstream. 


e. Go to step 6. 



5. 6 positive indicates expansion wave required. 

a. Compute the total pressure P tot based on upstream properties: 


tot 


P f"l M 2 

u L 2 u_ 


Y-l 




M U " q u' ' P 'u 


(3.21) 


where subscript u indicates values immediately upstream. 

b. Compute the corrected Mach number M c for deflection angle 6 and 
Mach number M 2 = q 2 / (YP/P) using Prandtl-Meyer Eq. (8.13) from 
Ref. (3-10): 


6 * 1 


jtan- 1 

(^)(M 2 -D 

- 4 


{tan* 1 ’ 
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(3.22) 


(Solution by Newton-Raphson iteration.) 


c. Compute the corrected static pressure P c using isentropic 
relationships: 


P = P 1 . 
c tot 2 c 


l+^M 2 




(3.23) 
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d • Compute the corrected static density p c using isentropic 
relationships: 


\ 


P 

c 



e. Go to step 6. 

6. Compute the corrected velocity magnitude q c using total enthalpy 
relationship: 


4 


2[h - <a> ( if> 1 


7. Compute the tangential components of velocity u>j, v^, and 

by discarding the normal components and calculate the tangential 
velocity magnitude qj: 


u = u - (q . n)' n 
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v T = v-(q.n)n 
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8. Compute the corrected velocity componenst u c , v c , and w c by 
scaling the tangential components by the ratio (q c /qx) : 



(3.27) 


The above procedure varies somewhat from that described by Abbett (Ref. 
3-8) but retains the essential features. Unlike Abbett's approach which 
applied the boundary condition once after the predictor step and then dis- 
pensed with the corrector step, the GIM code steady state marching algorithm 
applies the above procedure after both predictor and corrector steps. The 
boundary geometry of the current station is used on the predictor step since 
the corrected flowfield boundary values will be used with metric and flow 
field information from the current station to compute the corrector step. 
After the corrector step has advanced the solution to the next station, then 
boundary geometry from the new station is used to correct the flowfield 
boundary values at the new station. 

3.2 USER'S GUIDE 

3.2.1 General 

As detailed in previous sections, the GIM code steady state marching 
algorithm employs progressively assembled finite difference analogs during 
the integration process. This procedure differs considerably from that 
employed in the standard GEOM-MATRIX/MATQP-INTEG GIM code sequence. This 


3-16 



resulted in the development of two new G1M code modules, GEOMD and LNTKGII, 
for use in steady state marching problems. GEOMD is a specialized version 
of the standard GIM code GEOM module designed for use exclusively with 
integration modules that utilize progressively assembled finite difference 
analogs. INTEGH is the steady state hyperbolic integration module for the 
Euler equations. These two new modules are used together without any inter- 
vening matrix module. Both modules are described below together with a 
summary of the required input. Only those input items that differ from the 
standard GEOM and INTEG input are described in detail. For a detailed 
description of other input items and conventions, the user is referred to 
previously published GIM code literature (Refs. 3-11, 3-12, and 3-13). 

3.2.2 Dynamic Dime nsioner for GEOMD (DYNMAT) 

GEOMD uses the same dynamic dimensioning module as does GEOM, 
namely DYNMAT (Ref, 3-13). All input is identical. The user should, 
however, set up the dynamic dimensioning input data as if running a quasi- 
parabolic (QP) problem. 

3.2.3 Geometry Module (GEOMD) 

The geometry module used with the GIM hyperbolic Euler solver (GEOMD) 
is a modification of the standard GIM GEOM module. Due to the use of pro- 
gressively assembled difference analogs, all of the time consuming element 
matrix computations were removed, and no MATRIX module is required. The 
output has been made more readable and completely mirrors the input. 

The purpose of the GEOMD module is two-fold: (1) construct the com- 

putational mesh via trivariate blending, and (2) provide the integration 
module with connectivity information for each node (i.e., an array which for 
each node gives the node number^ of its nearest neighboring nodes along the 
three coordinate lines). 
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The output information is essentially the same as in other versions 
with two minor exceptions. First, the boundary node counter which pre- 
viously counted nodes with non-interior boundary flags (i.e., IB ss 9) now 
counts non-interior nodes and nodes with a connectivity which differs from 
the natural numbering scheme within any zone. This number can still be 
safely input as MNB in the dynamic dimensioner for the integration module. 
Secondly, these interior nodes with irregular connectivity are assigned a 
new boundary condition flag of IB--9, This is ignored in the integrator, 
but allows these nodes to be traced in the GEOMD output. 

GEOMD creates only a formatted FILE20 output file which contains 
geometry and connectivity information with several user input options. 
These are detailed in the following discussion of the input changes. 

Chart 3-1 gives an abbreviated description of the GEOMD input data. 
The input changes from previous documentation are given in the following. 

Card 2: IFMT replaces ISTEP. IFMT controls the 

output to FILE20 

IFMT^O Full GIM code output with coordinates, flow 
angles, normal components and boundary con- 
dition flags 

«1 Only Cartesian coordinates 

>1 Cartesian coordinates and boundary flags 

IMATRX«0 Connectivity array is computed and appended 

to FILE20 at each output record 

*1 No connectivity array. 
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Card Type 


1 HEADER(I)*, 1=1,72 

(12A6) 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 


NZONES, IDIM, IFMT, IMATRX, IMATE 
(515) 

IWRITE, LWRITE, NWRITE 
(315) 

NSECTS 

(15) 

'MAPE(I) , 1=1,12 
(1215) 

MAPS (I) , 1=1,6 
(615) 

(IBWL(I), 1=1,6), ITRAIN 
(715) 

(NNOD(I) , 1=1,3), (ISTRCH(I) , 1=1,3) 

(615) 

DIVPI(l) , 1=1,3 
(3E10.4) 

[AETA(J ,1) , 1=1, NNOD(J) ] , J=l, IDIM 
(8E10.4) 

[ (AC(I,K, J) , 1=1,8), K=1 , 5 ] , J=l,4, or 12 
(8E10.4) 

[AS(I, J) , 1=1,8], J=l,6 
(8E10.4) 

(PT(I, J) , 1=1,5), J=l,4, or 8 
(8E10.4) 

[ (PMAX(I,K, J) 1=1,5), ETAMAX(K, J) , K»l,4], 

J=l,4 or 12 
(6E10.4) 

Chart 3-1 - GEOMD Input Guide 


See GIM documentation for explanation of FORTRAN symbols, NASA CR 3157 and 
CR3369. 
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Ca rd 3 : 


Output Control 


IWRITE and LWRITE are now incremental 
indices. This aids in limiting connectivity 
and debug output. 

IWRITE=0 No intermediate print. 

=N Intermediate print for every N*-* 1 node 

LWRITE=0 No connectivity printout 

=N Print connectivity array for every N c ^ node 

Finally, notice that the card carrying the analog choice flags has been 
removed since the analog calculations are no longer performed in this module. 

The GEOMD module is FORTRAN 5 compatible and can be run on any serial 
computer. In particular, the CYBER 730-Z machine at NASA-Langley can be 
used for most two-dimensional problems and three-dimensional problems with 
no more than approximately 900 nodes in a cross plane. 

3.2.4 Dynamic Dimesioner for INTEGH ( DYNDIH ) 

INTEGH uses the same dynamic dimensioning module as does INTEG, namely 
DYNDIH (Ref. 3-13). All dynamic dimensioner input, which should be set up 
as if running a quasi-parabolic (QP) problem, is identical except for the 
last input item, IDYN. This variable must be set to 1 for use with INTEGH, 
as opposed to the value of 0 used with the standard INTEG module. 

3.2.5 Integration Module (INTEGH) 

The GIM code steady state hyperbolic integration module, INTEGH, is 
used to integrate the Euler equations by the methodology described pre- 
viously. The equations are integrated in space starting from the first 
upstream station (or plane) with a user-supplied flow field initialization. 
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The solution is then progressively "marched" downstream one plane at a time, 
the solution on each plane serving as the initial value surface for the 
integration to the next plane. The geometry and nodal connectivity inform- 
ation supplied by GEOMD via FILE20 is used to progressively assemble the 
finite difference analogs as required for each plane. 

Since the geometry information is constructed by GEOMD before the inte- 
gration process begins, it is quite possible that the spatial step size 
between the nodes of two planes may exceed that required for stability with 
the given flow field. The user can, therefore, direct INTEGH to monitor the 
spatial step size for stability considerations using the analysis of section 
3.1.2 (see METHOD on card H2 below). Then if INTEGH detects that the 
spatial step size between two primary planes (supplied by GEOMD) is too 
large for stable integration at one or more nodes, it will supply inter- 
mediate "ghost" or secondary planes by interpolating between the primary 
planes. This ensures that the integration can continue in a stable manner 
with the given geometry. If the user is certain that the geometry and 
anticipated flow field conditions are such that the spatial step size will 
always be within stable limits, then INTEGH can also be directed to cir- 
cumvent the stability analysis at each plane. However, limited experience 
has indicated that the calculations will fall apart almost immediately if 
the user is incorrect and the spatial step size is too large for stable 
integration. 

In some instances, the integration process may not be able to success- 
fully march to the next plane due to the presence of non-hyperbolic con- 
ditions, e.g., subsonic flow, flow reversal (separation), etc. INTEGH will 
normally detect these non-hyperbolic conditions, print an error message, and 
terminate the calculations. The user can, however, direct INTEGH to attempt 
to march through these regionq (see NHPLIM on card H2 below). This can 
sometimes be accomplished successfully if the region of non-hyperbolic con- 
ditions is not too large or too severe in content. 


c - < 
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Because It is a spatial marching integrator as opposed to the time- 
dependent integration of the standard INTEG module, the boundary condition 
treatment in 1NTEGH is of necessity somewhat different than that found in 
INTEG (see Section 3.1.4 for details). Because of this different treatment, 
all solid boundaries are treated as free-slip/tangency surfaces. Therefore, 
the user should note that INTEGH will only recognize and enforce boundary 
conditions of type 1 (axis node), type 3 (three-dimensional corner node), 
type 4 (free-slip/tangency node), and type 8 (free boundary node). Type 2 
boundary conditions (no-slip/stagnation node) are treated as type 4. Type 0 
boundary conditions (constant node) are not recognized at all and must not 
be used. In particular, the initial or inflow plane qnd the last or outflow 
plane should be treated with a type 9 boundary condition (no constraints) 
and not type 0, type 8, etc. 


The input data required for INTEGH are summarized in Chart 3-2 . In 
much the same manner as with a quasi-parabolic integration run, the user 
supplies INTEGH with the geometry information from GEOMD on FILE20, the 
initial flowfield values on the first plane or inflow surface, and other 
necessary control variables. Many of these input variables for INTEGH are 
identical to those required for standard INTEG runs and are described in 
detail in Refs. 3-11, 3-12, and 3-13. Therefore, only those input variables 
peculiar to INTEGH are described in detail below. 


Card H2: 

METHOD - 1 

- 2 


- 3 


METHOD controls the mode of operation of INTEGH. 

No check of spatial step size for stability. 
Cross planes must be orthogonal to the x-axis. 

Spatial step size checked for stability on each 
cross plane - secondary planes used if neces- 
sary. Cross planes must be orthogonal to the 
x-axis. 

No check of spatial step size for stability. 
Cross planes need not be orthogonal to the 
x-axis. 
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Card Type 


Parameter/Format 


HI 

H2 


H3 

H4 

H5 

H6 

H7 

H8 

H9 


ITITLE(I), 1=1,80 

(80A1) 

I DIM , METHOD , I TMAX , I TS AVE , I START , 10TYPE , IUNITS , 
1VISC.NHPLIM 

(915) 

DXFAC,EITOL 

(2E10.0) 

NPLT(I) ,MI2(I) ,MI3(1) ,NATT(I) ,L2(I) ,L3(1) , 
I=1,NZ0NES+1 

(615) 

GAM,WM,RK,EMU 

(4E10.0) 

KC ( I ) , 1=1,6 

(6A5) 

NJ , I NC , NTOT , ITAN , I TYPE 
(515) 

RI,UI,VI, WI,PI 

(5E10.0) 

N1 , IC , NT , NP1 , NPL , LSP 
(615) 


Chart 3.2 - INTEGH Input Guide 
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= 4 Spatial step size checked for stability on each 

cross plane - secondary planes used if neces- 
sary, Cross planes need not be orthogonal to 
the x-axis. 

ITMAX limits the number of iterations used to solve equations (3.5a, 
3.5b) when cross planes are not orthogonal to the x-axis. ITMAX is 
ignored if METHOD .LT. 3. 

ITMAX = Maximum number of iterations permitted. 

Nominal value « 5. 

ITSAVE determines whether or not the flow field will be saved on 
FILE22 . 

ITSAVE « 0 Flow field not saved on FILE22. 

* 1 Flow field saved on F1LE22. 

ISTART selects a cold start run or a restart run. 

ISTART = 0 Cold start run. 

= N Restart run beginning at plane number N. Flow 

field data read from FILE22. 


IVISC selects the use of pseudo-viscous terms to assist in the 
capturing of strong shocks. 

IVISC =0 No pseudo-viscous terms used. 

= 1 Pseudo-viscous terms used to assist capturing 

of strong shocks. 

(Note: This option which uses a variation of 
Harten’s method of Ref. 3-14 is not well 
exercised and should not be used at this time.) 

NHPLIMM selects the number of cross planes in which non-hyperbolic 
conditions can occur before terminating the calculations. 

NHPLIM = 0 Terminate calculations upon encountering 

non-hyperbolic conditions. 

~ N Up to N cross planes with non-hyperbolic con- 

ditions can be encountered before calculations 
will be terminated. 
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Card H3: 


DXFAC is a safety factor (analogous to DTFAC 
used in INTEG) applied to the spatial step size 
calculated from the stability analysis. DXFAC 
is ignored if METHOD = 1 or 3. Nominal value = 
0.9 - 1.0. 

EITOL is a relative tolerance used to terminate the iteration 
procedure used in equations (3.5a, 3.5b) on cross planes which are 
non-orthogonal to the x-axis. EITOL is ignored if METHOD - 1 or 2. 
Nominal value = 1.0E-05. 


Card H5: GAM is ratio of specific heats (same as GAMS1 

in INTEG) 

WM is the molecular weight of the fluid (same as WM1 in INTEG). WM 
can be left blank in which case RK must be supplied. 

EMU is a multiplicative coefficient used with the pseudo-viscous 
option (IVISC = 1). Nominal value = 0.5 - 1.0 

(Note: This option which uses a variation of Harten’s method of Ref. 
3-14 is not well exercised and should not be used at this time.) 


Card H9 : N1 is the node number of the first nodal point 

on this card to be printed during ouput. N1 = 
-1 terminates the reading of cards of type H9. 

IC is the nodal number increment to Nl at which output is desired. 

NT is the total number of nodes to be printed by this card type H9. 

NP1 is the first plane to be effected by the values of Nl, IC, and 
NT on this card. 

NPL is the last plane to be effected by the values of Nl, IC, and NT 
on this card. 

LSP selects the printing or suppression of secondary plane output 
associated with planes NP1 through NPL. LSP is ignored if METHOD = 

1 or 3. 

LSP = 0 

= 1 


No secondary plane output is printed. 

All secondary plane output associated with 
planes NP1 through NPL is printed. 
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3.3 SAMPLE CALCULATION: SHOCK EXPANSION TEST CASE 


Figure 3-2 illustrates the details of a sample calculation performed with 
the GIM code steady state spatial marching algorith (INTEGH). This test case 
consists of uniform Mach 2.4 flow in a two-dimensional channel which en- 
counters a 14.04 deg ramp followed by a complementing 14.04 degree turn back 
to the initial flow direction. An oblique shock wave with pressure ratio of 
approximately 2.3 forms at the intersection of the ramp and lower channel 
wall. This shock wave continues downstream and after reflecting off of the 
upper channel wall interacts with the expansion fan formed where the ramp 
turns back parallel to the upper channel wall. Figure 3-3 shows the computa- 
tional grid consisting of 85 cross planes with 41 modes per plane. Figure 3-4 
presents the velocity vectors computed by INTEGH. The dense packing of the 
velocity vectors highlights the well-defined shock wave. The expansion fan 
although present is less clearly discernable. Both the shock wave and the 
expansion fan can be seen, however, in Fig. 3-5 which is the same velocity 
vector plot as Fig. 3-4 but with every other vector masked out. Figures 3-6, 
3-7, and 3-8 show the pressure, Mach number, and density contours, respec- 
tively, for this shock-expansion test case. The pressure and density contours 

have been normalized by their inflow values. Computational time for this case 

-4 

was approximatey 10 CPU sec per node. Due to the short vector lengths 
encountered (41 nodes per cross plane), it is anticipated that INTEGH will be 
even more efficient for larger problems. 
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Fig. 3-2 - Shock Expansion Test Case (Diagram of Solution) 



Fig. 3-3 - Shock Expansion Test Case (Computational Grid) 
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Fig. 3-4 - Shock Expansion Test Case (Velocity Vectors - Dense) 



Fig. 3-5 - Shock Expansion Test Case (Velocity Vectors - Sparse) 
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Fig. 3-8 - Shock Expansion Test Case (Density Contours) 
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4. INVESTIGATION OF A SOLUTION-ADAPTIVE GRID ALGORITHM 
FOR THE GIM CODE 

4.1 INTRODUCTION 


The starting point of a numerical scheme for a digital computer is an 
array of finite points used to discretize the continuum. The treatment of 
fluid flows by these finite differencing methods has advanced tremendously 
in recent years due to production of advanced computing machines such as the 
CDC CYBER 203 and the CRAY. In a sense the machine advancements have pre- 
ceded the methods advancements. The pacing item appears to be the very 
starting point itself, the generation of a finite mesh of points on which to 
write the discrete equation set. 

We will not attempt to give a complete literature review in this re- 
port, as a separate report would be required for the extensive number of 
papers. The NASA-sponsored meetings at Langley (Reference 4-1) and Ames 
(Reference 4-2) contain articles of significance. The AlAA-sponsored 
meetings on Computational Fluid Dynamics (References 4-3 through 4-5) con- 
tain some articles on adaptive and multi-grid methods. The symposium, 
organized by Mississippi State University and sponsored by AFOSR and NASA 
(Ref. 4-6) provides a timely forum for this important topic of grid 
generation. 

The GIM (Ref. 4-3) code solves the multi-dimensional Navier-Stokes 
equations for arbitrary geometric domains. The geometry module in the GIM 
code generates two- and three-dimensional grids over specified flow regions, 
establishes boundary condition information and computes finite difference 
analogs for later use in the numerical integration solution module. The 
grid generation technique can be classified as an algebraic equation 
approach as opposed to a differential equation approach. The grid algorithm 
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uses multi-variate blending function interpolation of vector-valued 
functions which define the shape of the edges and surfaces bounding the flow 
domain. By employing blending functions which conform to the cardinality 
conditions, a unit cube can be mapped onto the flow domain, thus producing 
an intrinsic coordinate system for the region of interest. The intrinsic 
coordinate system facilitates grid spacing control to allow for nonuniform 
distribution of nodes in the flow domain. 

The current state of the art in Computational Fluid Dynamics relies on 
grids that are fixed in time, i.e., stationary Eulerian frames. Optimum use 
of computer storage and CP time cannot be made a priori on a fixed grid. 

The advanced algorithms will require time-dependent grids that adapt 
themselves to the physics of the flow. The grid movement terms are well 
known and a number of codes already exist with these terms programmed. What 
is missing is a methodology, even an algorithm, to cluster the points, i.e., 
supply values for these known terms in the equations and to control movement 
of the points. 

The objectives, then, of this research are: 

• To fill some of the information gaps which exist in adaptive 
grid generation technology 

• To provide a consistent treatment of grid generation and error 
control on physical domains, with nonuniform spacing 

• To develop a practical algorithm for a self-adapting grid for 
Navier-Stokes solutions. 


Our development will rely heavily on the General Interpolants Method 
(GIM) methodology, but the formulation and algorithm will be GXM independent 
for use by the general community. 
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4.2 FORMULATION 


This work begins by considering the Euler equations in two space dimen- 
sions and time. Extension to viscous flows and the Navier-Stokes equations 
is straightforward. The basic equation set in x,y coordinates is 


9U .. 3E . 3F 
3t + Sj + 97 


0 


u 



(4.1) 


r pu 1 


- pv 

2 



pu + P 


pvu 


F - 

2 , „ 

pu V 


pv + P 

(p£ + P) U_ 


_ (p£ + P) v_ 


Here x,y,t are the spatial coordinates and t is time. The flow variables 
are p * density, u,v = velocity components in x,y directions, respectively, 
£= total energy and P = static pressure. 

A set of "boundary conforming" coordinates is generated by a geometry 
module GIM code to simplify the boundary condition procedures and difference 
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scheme. A computation domain is constructed with uniformly spaced rectan- 
gular coordinates and the physical coordinates are mapped onto this domain. 


(x,y,t) 


(£» n» t) 



Domain 



Computational 

Domain 


Insertion of a general curvilinear coordinate transformation, keeping all 
time-dependent terms, results in the following equation set: 



* * 
3E_ 3F_ 

H 


= 0 


where 


if = JU 

E* = y n [E - X T U] - X n [F - y T U] 
F* = x^[F - y T U] - y ? [E - x^] 


( 4 . 2 ) 
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Here J is the determinate of the transformation Jacobian 


J 


x y - 

rn 


x y r 
n K 


( 4 . 3 ) 


and x_ , y_, x , y are the metric coefficients and 

K. K n n 

V 

x t 57 • \ ‘ 

are the coordinate movement velocities in the x,y directions, respectively. 

Thus the coordinates in physical space can move in time according to 

the prescription for x , y . This can be from a moving boundary or from 

T T 

interior point movement via an adaptive grid algorithm. At each time step, 
the physical grid is moved and remapped onto the fixed computational grid 
via the metrics and the Jacobian. The finite difference program for solving 
the equations is thus unaltered, and the effects of coordinate movement is 
taken into account in the metric data. 

Generation of these metric data, and in particular the Jacobian J, must 
be done with extreme caution. This is discussed in detail by Thomas and 
Lombard (Reference 4-7) where they introduce the Geometric Conservation Law 
(GCL). In this procedure, a differential equation is solved for the 
Jacobian J using the same difference scheme as used on the flow equations. 
This results in total consistency of numerics and preserves uniform flow 
exactly. In an adaptive grid scheme this will be very important to a good 
working code. We thus append the following equation to our set. 


h. 


To complete the formulation, we need a prescription for determining the 
coordinate movement x T , y T . For the first version of the GIM adaptive grid 
code, we use the scheme of Rai and Anderson (References 4-8 and 4-9). In 
one space dimension, this scheme is simply written as 



(4.5) 


where e^ is a measure of the error at node j, 
nodes i and j, and K,M are chosen constants, 
scheme in one of three ways: 


r^j is a distance between 
The current code uses this 


• Coordinate movement in x direction only 

• Coordinate movement in y direction only, and 

• Coupled movement in both directions. 

The functions ej present one of the major difficulties in using this 
scheme. Accurate determination of the truncation error in our scheme has not 
been done yet. In the current code, the functions are set to the local 
pressure gradient, 


e 


j 


9P 

at 


i 


Control of the grid speeds, , y^ is also a concern with this 
scheme. Control of the speed is currently done with the constant K, chosen 
a priori in the manner of Rai and Anderson (References 4-5 and 4-6). 
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The MacCormack explicit , unsplit , finite difference scheme (Reference 
4-10) is used to solve all differential equations. Boundary conditions con- 
sist of fixed in-flow, supersonic one-sided out-flow and solid wall 
free-slip (tangency). The grid velocities are forced to obey wall tangency 
so that the boundary points stay in the domain. 


The computation proceeds as follows: 

• An initial grid is prescribed on physical space by the grid 
coordinate generator. 

• These coordinates are mapped (numerically) onto the computa- 
tional domain and the metric data calculated. The Jacobian J 
is initialized by the algebraic definition (Eq. (4.3)). 

• The grid velocities are set to zero. 

• Equation set (4.2) is formed and solved for U* by the 
MacCormack method. Note that U* * (pJ, puJ, pvJ, p£j,j). 

• Boundary conditions are applied. 

• The vector U* is decoded for p, u, v,£, J. 

• The pressure P is determined from the ideal gas law. 

• Grid velocities are updated from Eq. (4.5). 

• The coordinates x,y are moved using a first-order scheme 

x n+1 * x n + x n Ax, 
t ’ 

yH +1 * yn 4 . y n 

• New metrics are determined, x , x , y. , y from the known 
coordinate transformation. s H c, H 


• Equations (1) are formed (now with x^_ / 0, y^ / 0). 

• The solution process is repeated until convergence is achieved 
to steady state or for a fixed number of time steps. 
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4.3 RESULTS 


An example case, simple in nature, is shown here to illustrate that 
this formulation does work. This example case is a standing shock wave in a 
duct . 


In Flow 



Out Flow 


The shock was placed in the axial center of a fixed grid of 26 x 5 points 
(5 points in y direction for simplicity). The code was run until the grid 
velocity x^ dropped six orders of magnitude from its peak value. 

Converged results are shown in Figure 4-1. The left figure is the 
MacCormack solution for a fixed grid, i.e., x T * 0 for all T. The right 
figure shows the result for an adaptive grid. The plots are pressure versus 
x with the grid point positions shown at the bottom for reference. As seen 
here, the grid motion resulted in points being moved closer to the shock 
discontinuity and also resulted in a smoother calculated shock, i.e., the 
overshoot is almost gone. 

Figure 4-2 shows the same problem run with 51 axial points to test the 
algorithm on a finer mesh. The left figure shows virtually the same results 
for grid movement, but of course, a steeper calculated shock. The right 
plot on Figure 4-2 is an example of this problem when the constraint on x 
was relaxed. The grid speeds get too fast, the points cross one another, and 
the entire solution goes unstable. 
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Pressure vs X Pressure vs X 



Fixed Grid Adaptive Grid 

^ Figure 4-1. One-dimensional shock tube problem-adaptive grid solution 

i 



Converged Caee Unstable Case 

(Fine Mesh) (Too Little Constraint on x T > 


Figure 4-2. One-dimensional shock tube problem-adaptive grid solution 
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4.4 EXTENSION TO HIGHER DIMENSIONS 


The work on adaptive grids done under the NASA-Langley contract was 
restricted to the one-dimensional development. This work is currently being 
extended under Lockheed's Independent Research Program. Two- and 
three-dimensional formulations have been accomplished and coding is now in 
progress. Documentation of the multi-dimension work will be done under the 
Independent Research Program. 
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5. TURBULENCE MODELING IN THE GIM CODE 

5.1 INTRODUCTION 

This section updates the modifications and additions to the GIM code 
turbulence models since their last documentation (Ref. 5-1). The algebraic 
eddy viscosity model of Baldwin and Lomax (Ref. 5-2) has been extended to 
two-dimensional flows bounded by more than one wall. The option is also 
provided to apply the model selectively in different regions of the flow. 
These changes are discussed in the next section. The third section details 
the final form of the wall function boundary conditions chosen for the tur- 
bulent kinetic energy (TKE) differential equation model, and describes their 
use and initialization. 

The turbulence models exist as files of UPDATE directives to the basic 
GIM code INTEG module. This approach was adopted so as to not adorn the 
basic module with often unused subroutines and so that the models could be 
coded in a manner compatible with the INTEG module but not necessarily with 
each other since only one is used at a time. 

5.2 IMPROVEMENTS TO THE BALDWIN/ LOMAX MODEL 

The Baldwin/Lomax algebraic eddy viscosity model in the GIM code has 
been extended to allow the user to apply the model within as many as ten 
different regions along a cross-flow line of nodes. This approach allows 
the user to compute in regions with more than one wall and/or selectively 
apply the model to different regions of the flow domain. As previously 
documented, the model can be used for axisymmetric and two-dimensional flows 
in either the quasi-parabolic or elliptic modes of GIM calculation. It is 
assumed that one set of grid lines lies approximately in a cross-flow 
direction (U^ for QP calculations). 


1 
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These features are implemented by modifying DATA statements within the 
UPDATE file itself. The variables which need to be changes are 

JREG * The number of regions along a line of nodes 

J1RG(I) - The starting node for Region I along a line 
of node 8 

JLRG(I) « The number of nodes included in Region I. 

These variables are presently set for one region extending over the full 
line of nodes. J1RG is a local node number; i.e. , 1 <_ J1RG <_ number of 
nodes in the cross-flow n direction. 

Experience indicates that the model must be modified to some extent for 
almost every calculation performed. This is a result of the coupling of the 
generality of the G1M methodology and the specific nature of most algebraic 
turbulence models. To provide further versatility (especially in geometric- 
ally complex problems requiring more than one zone); the user can modify the 
defining statements of the following two variables: 

JFRST = The global (GEOM-assigned) starting node number in 
the region along the line of nodes being computed. 

NNY ■ The number of nodes to be included in the computation. 

In the general case NNY ■ JLRG. 

A simple example of flow in a region between two walls may serve to 
clarify the use of these regions. The problem is illustrated in Fig. 5-1. 

is the streamwise direction and is the cross-flow direction. Here 
we choose to divide the flow into two regions at the 25th node along the 50 
node u 2 line. The modified DATA statements in the update file for this 
case should appear as 
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Node JFRST+49 



Fig. 5-1 - Flow in a Two-Dimensional Region Bounded by Two Walls 


DATA JREG/2/ 

DATA J1RG/1,26,8*0/ 

DATA JLRG/25,25,8*0/. 

In the event that the minimum velocity within the region, which the 
model uses as a reference location, does not occur at one of the end points 
of the region (for example, this situation could occur in wake flows) the 
model will subdivide the region at the minimum point automatically. 

Finally, a word of caution. Although the above may seem complex, it 
has been found to be a coding arrangement ammenable to the modifications 
which seem necessary with almost every new problem. It was coded in this 
fora so as not to disturb the input to the baseline model and at the same 
time to require considerable usqp: preparation prior to modification. 
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5.3 BOUNDARY CONDITIONS FOR THE GIM/TKE MODEL 


The boundary conditions for the two-equation turbulent kinetic energy 
TKE model in the GIM code have been finalized. They consist of wall func- 
tion descriptions of the behavior of the momentum, turbulent kinetic energy, 
k, and the rate of dissipation of the turbulent kinetic energy, e. These 
conditions override the type 2 (no-slip) boundary conditions in the code 
when the TKE model updates are applied. The model can be used for axi- 
symmetric, two- and three-dimensional elliptic or quasi-parabolic calcula- 
tions. This section will describe the form of these boundary conditions and 
their application in the code. To date, no test cases have been run using 
this model and the need for artificial damping to stabilize the time- 
iterative solution of these equations has not been studied. The form of the 
damping terms has not been derived. 

Wall function boundary conditions for differential equation turbulence 
models allow these more complex and storage-intensive models to be computed 
without resolving the behavior of the laminar sublayer which necessitates 
very fine meshes and large number of nodes. It is assumed that the first 
node off the wall in such calculations lies either in the laminar sublayer 
or the logarithmic "law-of-the-wall" region. The wall functions in the 
GIM/TKE model are modeled after those of Taylor, Thomas, and Morgan (Ref. 

5-3) written in a compressible form. These functions are applied at type 2 
boundaries for the momentum, k and £ equations. The wall functions are 
given in Fig. 5-2. For the momentum equations the computed flow direction 
is constrained to be tangent to the wall and the magnitude of the momentum 
is computed from the wall function. 

The variable viscosity flags are input as follows for the GIM/TKE model: 
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Momentum: 

+ 

pn ; 

o < n + < 5 

+ 

pq 

= | P(5.0 £nn + - 3.05) (t w /|t w |); 

5 £ n + < 30 


P (2 . 5 £nn + + 5.5) (x / [ T |) ; 
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Dissipation of TKE 

(pe) w - C y (k w ) 3/2 /< n 


where 



T 

W 

u L (n . V) (t . v) n + - (pn/y L ) y|i w | /p 

+ 

q 

pq/V|x w |/p q - 1^1 



K - 0.42 



Fig. 5-2 - GIM/TKE Wall Functions 
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0.09 


VF(3) = C y = 

VF(4) = C £l - 1.44 

VF(5) = C e2 - 1.92 

VF(6) = Pr t = 0.9 

VF(7) « y x /y L “ 10 to 100 

VF(8) » q'/q * a few percent 

where C y , and are constants in the model, Pr t is the turbulent Prandtl 
number, y^/y^ Is the ratio of turbulent to laminar viscosity, and q'/q is the 
ratio of turbulent fluctuations of velocity to the magnitude of velocity. 

The latter two quantities are used in initialization of the turbulence quan- 
tities. Throughout the flow k and £ are initialized from input flow quan- 
tities as 


k - j [q . (q’/q)] 2 

£ * P C y k 2 /y x 


where kp is determined from y^ and VF(7). These quantities are insuf- 
ficient to initialize the magnitude of the momentum at the wall points. 

Here we must assume the normal distance of the first node point from the 
wall. In the code it is assumed that the first node is one-tenth the normal 
grid spacing from the wall, i.e., 


n w " (n x Ax + n y A y + n z Az )/ 10 


where n^, n^, and n z are the components of the normal to the surface 
at the wall node and Ax, Ay and Az are the grid spacings at that node. The 
dimensionless speed becomes 
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from this the dimensionless normal distance from the wall can be calculated 
from the wall functions as 


+ 

n 


f q + 


exp[ (q + 3.05J/5.0] ; 


exp[(q - 5 .5) /2 .5] 


0 < q + < 5 
5 ^ q + < 14 
q + > 14 . 


Finally, the dimensional speed at the wall becomes 


q 


+ , 

Um H + 

m T 1 q 

p (y T /y L ) 


The energy Is Initialized consistant with the input pressure and density and 
the above speed. It is necessary for the user to input type 0 (fixed) 
boundaries in a consistant fashion. 
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GENERATION OF RATE CONSTANTS FOR A GLOBAL 
HYDROGEN- AIR COMBUSTION MODEL 


Appendix A 


A.l INTRODUCTION 

The analysis of non-reacting flow fields is still a difficult task, 
even with modern computational techniques. The task is even more difficult 
when reacting flows have to be considered. Realistic modeling of finite 
rate reaction systems usually requires the consideration of numerous chem- 
ical species, and often even more numerous reaction processes describing the 
interaction of the various species. This, in some cases, drastically in- 
creases computational time, and always increases program storage require- 
ments. It is therefore of interest to devise procedures which can be used 
to establish reduced or global reaction models which are capable to faith- 
fully reproduce the results of detailed and complete models, in spite of 
using only a limited number of species and reaction equations. Provided 
that such a global mechanism can be defined, which is probably more a matter 
of proper judgement than of computation, a systematic procedure to evaluate 
the necessary global rate constants is described here. The specific case of 
hydrogen-air combustion in a supersonic flow is considered, but the pro- 
cedure should be applicable to arbitrary reaction systems. 

A. 2 METHODOLOGY 


The application of the methodology to be described assumes the avail- 
ability of a time-dependent finite rate reaction computer code which can 
compute the evolution of the species composition and the temperature as a 
function of time during a combustion process. A one-dimensional, inviscid 
pre-mixed calculation suffices for this purpose. 


A-l 


Let us assume that the full set of species and reactions to be reduced 

(or "globalized") is given by m = 1 NR reactions involving i = 1,...,NS 

species, i.e., 


yv a. lz Vv; 

^ i,m 1 "• — 2L j i,m 


(A.l) 


where the A. denotes the species chemical symbol and the v, and v! 

i r i,m i,m 

are the stoichiometric coefficients of the reactants and the products for 

the m*"* 1 reaction, respectively, k. denotes the forward rate constant 

th 1 ,m 

for the m reaction. 


The net rate of production for any species A^ participating in re- 
action m (i.e., v! - v. ^ 0) is then given by 

l , m l ,m 


w 




(Vl - V ) w 
i,m i,m m 


(A. 2) 


where 


• r i v' n 

w - k, I IT n c. 1 --] 

m f,m Li i K i i J 

* c,m J 


(A. 3) 


and where we have made use of the relation relating the forward and backward 
rate constants through the concentration equilibrium constant, viz. 


K 


c,m 


f ,m 

t, 

b,m 


(A.4) 


The total net rate of production for species from the set of m reactions is 
then given by i 
actions, i.e., 


then given by summing the contributions w. from (A. 2) over all m re- 

1 


T », 

Zj i,m 
m 


(A. 5) 
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Now let us assume that the above described full set (1=1,..., NS, m=l,.,.,NR) 

/\ 

Is to be replaced by an equivalent, but reduced set of j=l, ... ,NS species 

A * A 

and n=l,...,NR reactions, where NS < NS and, preferrably, NR« NR. The 
question that arises is, what are the rate constants for the postulated re- 
duced set of reactions? Assuming that by some reasoning a reduced set of 
species and reactions has been defined, the latter is described by 


£v< A . — £ V' A. 
j»n j *-£ j,n j 


and, as before, 


w. = (V - v. ) w 
J»n j,n j ,n n 


where 


w = k. 
n f 


r nc > n --i- n C>“1 

> n l-i j K c>n j j J 


(A.6) 


(A.7) 


(A.8) 


Again, the total net rate of production for species j from the reduced set 
of reactions is obtained as before by summing the contributions of all n 
reactions such that 


w 


j 


where AVj, n is defined as 






w 

n 




V 


j*n 


(A. 9) 


(A. 10) 


For j species comprising the reduced or "globalized" set, Eq. (A. 9) 
represents a set of j simultaneous linear equations for the unknown w n , 
i.e. , 




w , 


(A. 11) 
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where (Av ) represents a coefficient matrix, and w and w, represent 
j,n n j 

column vectors. 


In general, provided that n=j and that the coefficient matrix is non- 
singular the above set can be solved, either by successive elimination (tri- 
angulation) or by using Cramer* s rule, in which case the solution is 


w 

n 


(Av 


J»n 


) 1 w, 


(A. 12) 


where 




) 


j Av 



(A. 13) 


Since the reduced or globalized set is supposed to be equivalent to the 
original full set, we must prescribe that the net rates of production from 
the reduced set be equal to those of the full set for the species retained, 
i.e. , 


Wj - ( j*i=l , . . . ,NS) (A. 14) 

and, similarly, that the concentrations be equal, i.e., 

Cj - C ± ( j=i=l, . . . ,NS) (A. 15) 

Using these equivalency conditions, the left hand side of Eq. (A. 8) is 
determined by Eq. (A. 12), while the bracketed term on the right hand side of 
Eq. (A. 8) is determined by the known concentrations from the full set cal- 
culations. Hence we can solve for the unknown rate constant k, , i.e., 

f ,n 


k 


f,n 


w 


n 


v 4 


n - - 

j J K 


V 


n c. 


j.n 


c,n 


(A. 16) 
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where, according to Eq. (A. 15), the are the concentrations as deter- 
mined from the full set. 

The procedure is completed by plotting and curve fitting the k^ n as 
a function of temperature. Other parameters such as fuel to oxidizer ratio 
may have to be incorporated to generalize the applicability of the rate con- 
stants thus derived. 

Several problems can arise with the solution of the set of linear 
simultaneous equations represented by Eq. (A. 11), depending on the number of 
species j, and the number of reactions n, comprising the "globalized" set. 

If j > n, obviously the set is over-determined, and we have a choice of 
species production rates to be used to determine the reaction rates w^. 

If j < n, there are more unknowns than equations. The solution to this 
problem is to specify n-j rate constants, leaving j equivalent rate con- 
stants to be determined. Even if j*n, the coefficient matrix (v. ) is 

J» n 

not necessarily nonsingular. The reason for this is as follows. If we con- 
sider j species composed of i elements, we can generally write i algebraic 
equations for elemental conservation. This means that j-i rate equations 
are sufficient to compute the species composition of the gas, that is, only 
j-i equations are linearly independent. Considering the present application 
this means that again in Eq. (A. 11) n < j, which is the case discussed first. 

The specific example of hydrogen-air combustion described in the next 
section will clarify the application of the methodology just described. 

A. 3 HYDROGEN- AIR COMBUSTION 


The configuration chosen to test the general methodology was a rela- 
tively simple one, namely premixed, one-dimensional supersonic flow of a 
hydrogen-air mixture in a constant area channel. Detailed combustion 
calculations were performed using the ALFA code (Ref. A-l). Species 
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While 


considered consisted of N 2 > 0 2 > H 2 > H 2 O, OH, 0, H, H 2 0 2 and HO 2 . 

N, was considered to be inert concentrations of the remaining species were calcu- 
lated using a set of reactions for hydrogen-oxygen combustion adapted 
from those given by Beach, Mackley, Rogers, and Chinitz (Ref. A-2). The full 
set as used is given in Table A-l. 

On the basis of arguments presented by Rogers and Chinitz (Ref. A-3), 
the global model was selected to consist of two reactions, namely 

kf . 

h 2 + ° 2 2 OH (A. 17) 

kf 2 

2 H„0 2 OH + H 2 (A. 18) 

The second reaction actually being the reverse of the one given in Ref. A-3 
so to conform to the reaction rate input format of the ALFA code (since the 
reactions can proceed either way, the form chosen is immaterial). Species 
involved in the global model thus consist of N 2 , 0 2 , H 2 , OH, and 
H 2 0, i.e., five species as compared to nine in the full set, with a 
reduction in reactions from 18 to 2. 

Application of the previously described methodology to determine the 
unknown rate constants k^ ^ and k^ 2 is straightforward. 

From Eqs. (A. 8) and (A. 9) we find that 


w. 


‘£,1 < C H 2 % ' k ] (1 C 0H ) 


(A. 19 ) 


W 2 " k f,2 (C H,0 “ K , C 0H C H 9 J 
2 c,2 2 


(A. 20) 


See Appendix B, Table B-la and associated text for discussion. 
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Table A-l - DETAILED HYDROGEN-AIR COMBUSTION KINETICS MECHANISM 
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ORIGINAL page is 

0F POOR Quality 


and 




- "l + ”2 


w «* — 2 w* 

H 2 0 2 


w 


OH 


2 w 1 + 2 w 2 


(A.21) 


Since the net production rates of 0 2 , H 2 , H 2 o and OH are known from 
the full set calculations, we can determine w^ and w 2> As seen from Eq. 
(A.21) we have four equations for two unknowns, two of the equations pro- 
viding our unknowns directly, i.e.. 



(A, 22) 


Alternatively we can use the second and the fourth equation to obtain 


1 * 1 * 

u a — u — — tj 

1 4 OH 2 H 2 

1 * . 1 * 

W 2 = 4 W OH + 2 W H 2 

Thus, from Eqs. (A. 19), (A. 20), and (A.23) 


(A. 23) 



1 ’ 1 * 

4 W 0H ” 2 W H 2 



(A. 24) 
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(A. 25) 


4 U 0H + 2 U H. 

k e ■ ■ . . - - — ■ 

f 2 2 12 

1 n*’ _ p*' p 

H„0 K „ OH H_ 

2 c,2 2 

where the right hand sides are completely determined from the full set cal- 
culations. 

A. 4 RESULTS 

Results from the full set calculations are shown in Figs. A-l and A-2. 
Initial conditions were T = 1000K, p = 760 torr, M * 2.24 and an equivalence 
ratio of <|> = 0.25 (initial calculations with higher equivalence ratios pro- 
duced too large a heat release with choking of flow as a result). As seen 
from Fig. A-l, this case is characterized by a definite ignition delay (0 < 
x < 22 cm), before the temperature gradient shows any change (at constant 
velocity the distance of 22 cm corresponds to a delay time of 0.15 msec). 

The combustion front is centered around x = 29 cm, and combustion is basic- 
ally complete at x = 40 cm. 

Figure A-2 shows the species distributions through the combustion 
front. The two species not shown (H^ and H0 2 ) never reach relative 
concentrations of 10 \ Note that in the evaluation of k^ ^ and k f 2 
(the global rate constants) we ignore the existence of the minor species 
H 2 0 2 and H0 2 , as well as 0 and H, Including their transient peaks. 

Global rate constants k, , and k. » were evaluated simultaneously 
with the results just shown. They are plotted in Figs. A-3 and A-4, re- 
spectively. Note that k^ ^ is plotted versus T « T-T (initial). From 
the results it can be seen that k^ ^ increases by almost three orders of 
magnitude while the temperature increases by only one degree during the 
ignition delay period. It was found Impossible to fit this behavior with 
any of the standard rate constant formats (k * A.T exp(B/T)). As shown 
in Fig. A-3, k, ., however, displays a reasonable behavior when plotted 

i , 1 
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Fig. A-l - Supersonic Channel Flow with H^-Air Combustion (<Ji * 0.25) 
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Fig. A- 2 - Supersonic Channel Flow with H 2 ~Air Combustion (<fr » 0.25) 
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versus AT. plotted versus inverse temperature, is shown in Fig. 

A-4. The almost perfect straight line immediately suggests a fit of the 
form k ■ A exp(-B/T). Trial curve fits and calculations with the global 
model revealed that excellent results can be obtained from the global model 
provided that the global rate constants are curve fitted accurately in the 
low temperature regime (i.e., for the conditions during the ignition delay 
period). A close fit for the high temperature region is much less important 
as in this region the rate of production is essentially determined by the 
concentrations of the reactants being depleted. 

3 

The final best curve fits obtained are (in cm /part sec) 

k f 1 = 6.925 . 10 -17 AT 0,979 (A.26) 

k f 2 - 3.80 . 10" 6 exp(-116,000/RT) (A.27) 

with emphasis on accuracy in the low temperature regime. The validity of 
these particular rate constants is, of course, at least at this time, re- 
stricted to the conditions used in these calculations. 

Using the ALFA code, combustion calculations were then repeated using 
identical initial conditions but with the global model consisting of only 
five species (N 2> 0 2 , H 2 , OH, and H 2 0) and the two reactions (17, 18) 
with rate constants as given in Eqs. (A.26) and (A.27). The results are 
shown in Figs. A-5 and A-6 in comparison with the results from the detailed 
reaction model. As seen in Fig. A-5, results for the distribution of veloc- 
ity, temperature and pressure are in excellent agreement with the detailed 
results, and the ignition delay is modeled exactly, for all practical pur- 
poses. Final temperature (and pressure) is slightly higher than that from 
the detailed model, a result expected based on the discussion given in Ref. 
A-3. Species concentration distribution through the combustion front can be 
compared in Fig. A-6. The agreement is very satisfactory. 
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NS = 9, NR - 18 (Detailed Model) 
NS = 5, NR = 2 (Global Model) 



Fig. A-5 - Supersonic Channel Flow with H 9 -Air Combustion (<(> 






A. 5 CONCLUSIONS AND RECOMMENDATIONS 


A general methodology has been presented to systematically reduce de- 
tailed reaction models, employing numerous species and reaction rates, to 
global systems with a minimum number of species and reaction mechanisms. 
Global rate constants can be evaluated directly thus avoiding possible 
tedious trial and error procedures. 

The methodology was applied to the hydrogen-air combustion kinetics 
systems. Global rate constants were evaluated, curve fitted, and used to 
almost exactly reproduce the results obtained by using the detailed kinetics 
mechanism. The only real problem encountered was that of curve fitting one 
of the two global rate constants. It was found that, in order to model the 
ignition delay correctly, care must be taken to obtain an accurate repre- 
sentation of the global rate constants as a function of temperature in the 
ignition delay period. 

The global rate constants which were evaluated are obviously valid for 
the sample case presented, but not for other conditions (temperature, equiv- 
alence ratio, pressure, etc.). Further calculations should be performed to 
establish global (or compound) rate constants of more general applicability. 
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SUPERSONIC AIR-METHANE-SILANE/H 2 IGNITION* 


This work was performed in the period extending from 27 May 1981 to a cut- 
off date of 31 December 1981. 
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ABSTRACT 

Detailed kinetics mechanisms for silane-hydrogen and methane oxidation 
have been combined and adapted for use in computation of the characteristics 
of supersonic air-methane-silane/^ ignition and heat release. After 
initial mechanistic refinement interfacing with the presently meager data 
(1981), predictions made with Lockheed's ALFA code using the refined mecha- 
nism as presently constituted - in a constant pressure, one-dimensional flow 
regime - are at least within a factor of 2 of early shock tube measurements 
of ignition delay at Langley Research Center over the range 800 to 1000K for 
sllane-H2-02/N2 mixtures. Initial results of a parametric examination 
of theoretical ignition and reaction times in silane-CH^-air mixtures have 
been obtained. Use of silane as an ignition aid materially reduces ignition 
delay and combustion times for methane in air. 

INTRODUCTION 

The Inherent promise of silane as an Ignition aid In testing hydrogen- 
fueled, airframe-integrated, supersonic combustion ramjet engines was first 
reported by Beach et al. in Ref. B-l. Initial experiments and analytical 
research at NASA-Langley relating to the use of silane in scramjet ground 
testing is presented by Beach et al. along with an early SiH^/l^ combus- 
tion mechanism formulated to permit kinetic computations of hydrogen-air 
flows employing silane ignition aid. 

Use of hydrocarbons as fuel in potential scramjet applications is also 
of more than conceptual’ interest. As enumerated in the review of potential 
scramjet aircraft technology issues by Jones and Huber (Ref. B-2), 
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hydrocarbons suffer by comparison to hydrogen as scramjet fuel due to their 
lower heats of combustion per unit mass and longer ignition delays. Other 
systems aspects arising from their intrinsically higher density and ease of 
handling and storage as liquids combine, however, to make the choice more 
competitive. Future systems development may very well hinge on the crucial 
issue of fuel choice between hydrogen and hydrocarbons such as methane, 
methanol or conventional JP-type aviation fuels. The outcome of such compe- 
tition will clearly be effected by the utility of ignition aids such as 
silane in holding ignition and reaction times to manageable levels for these 
various fuels in flight regimes in which autoignition is not achievable. 

The present study has been undertaken with the goal of extending the 
work of Ref. B-l from the SiH^/I^ combustion system to the SiH^CH^ 
combustion system as a first step in the analytical exploration of the 
utility of silane as an ignition aid for testing hydrocarbon-fueled, 
airframe-integrated, supersonic combustion ramjet engines. Parametric com- 
putations were carried out to guide the experimental studies which seem 
likely to follow the favorable results of the present effort and initial, 
unpublished engine tests with hydrocarbons at GASL employing silane as igni- 
tion aid. During the course of the present study, results of initial shock 
tube measurements of ignition delays in SiH^/C^ mixtures diluted with 86 
vol.% N2 became available and it was possible to perform some mechanistic 
refinement based on these inputs, with the result that the SiH^/^ mech- 
anism of Ref. B-l has also been substantially updated through necessity to 
improve agreement with these data. 

REACTION MECHANISM 


A workable methane/hydrogen-silane combustion mechanism has been formu- 
lated. As presently constituted it is as shown in Table B-l, consisting of 
(A) hydrogen; (B) methane; (C) silane; and (D) methane-silane interactive 
oxidation mechanisms, discussed in brief below. Lockheed's ALFA computer 
code (Ref. B-3) has been utilized in this study for all kinetics computa- 
tions, in a constant pressure, premixed flow mode option. ALFA is a 
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Table B-l 

REACTION MECHANISM 


A. Hydrogen Oxidation Mechanism 

3 1 

Reaction (Reversible) Rate Coefficient, cm particle - Reference 

sec - "* units 


1. 

H 2 + ° 2 — OH + OH 

1.7 x 10" 11 exp-43, 000/RT 

1,7 

2. 

H + O z -*OH + O 

- 13 

1.0 x 10 T exp- 14,800/RT 

8, Text 

3 . 

OH + H., — *■ H z O + H 

1.8 x 10" 15 T 1 ' 3 exp-3, 650/RT 

8 

4 . 

O + H 2 —OH + H 

3.0 x 10“ 14 T exp-8, 900/RT 

8 

5. 

OH + OH — *> H-,0 + O 

1.0 x 10" 16 T 1 * 3 

8 


M 

, , ,„-26 _-2 


6. 

H + OH + M — *H 2 0 + M 

6.1 x 10 T 

9 

7. 

H + H + M — H 2 + M 

1.8 x 10" 30 T" 1 

1.7 

8. 

H + 0 2 + M — H0 2 + M 

5.8 x 10" 30 T" 1 

1,7 

9 . 

H0 2 + OH — H 2 0 + 0 2 

8.3 x 10" 11 exp- 1, 000/RT 

1,7 

10. 

H0 2 + H — H 2 + 0 2 

4.2 x 10“ 11 exp-700/RT 

1,7 

11. 

H0 2 + H —OH + OH 

3.3 x 10“ 10 exp- 1,800/RT 

1,7 

12. 

H0 2 + O —OH + O z 

8.3 x 10“ 11 exp- 1, 000/RT 

1,7 

13. 

ho 2 + ho 2 — h 2 o 2 + o 2 

3.3 x 10“ 12 

1,7 

14. 

HO Z + H 2“* H 2°2 + H 

5.0 x 10" 13 exp-18,700/RT 

1,7 

15. 

H 2 0 2 + OH — HC> 2 + H z O 

1.7 x 10 -11 exp- J, 900/RT 

10 

16. 

H 2°2 + H — OH + H z O 

8.3 x 10“ 10 exp-10, 000/RT 

10 

17. 

H 2°2 + ° ~~ OH + h °2 

3.3 x 10 -11 exp-5, 900/RT 

10 

18. 

H 2 0 2 + M —OH + OH + M 

2.0 x 10" 7 exp-45, 500/RT 

1,7 
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Table B-l (Continued) 


B. Methane Oxidation Mechanism 


Reaction (Reversible) 

(To be used with H, Mechanism) 

c -i _ 1 

Rate Coefficient, cm J particle 
sec" 1 units 

Reference 

1. CH. + M — -CH, + H + M 
4 3 

4.0 x 10 7 T" 3,7 exp- 103,200/RT 

11, Text 

2. CH. + O, — *CH, OH + OH 
4 2 2 

1.0 x 10“ 12 exp-43,000/RT 

12, Text 

3. CH + OH — -CH, + H,0 
4 3 2 

3.2 x 10" 18 T 2,1 exp-2,450/RT 

13 

4. CH 4 + H ^CH 3 + H 2 

3.7 x 10 T exp-8, 700/RT 

14, 15 

5. CH. + O — ►CH- + OH 
4 3 

1.1 x 10" 18 T 2 - 5 exp-7, 400/RT 

16, 17, Text 

6. CH . + HO, — CH, + H,0, 
4 2 3 2 2 

3.3 x 10 H exp-18, 000/RT 

18 

7. CH 3 + 0 2 — CH 2 OH + O 

1.2 x 10' 11 exp-25, 600/RT 

19 

8. CH 3 + 0 2 — CH 2 OH + O 

2.0 x 10“ 8 T~ 3 exp 1,000/RT 

31, Text 

9. CH 3 + OH — ~CH 2 OH + H 

7.0 x 10" 12 

19,20 

io. ch 3 + o — ch 2 O + h 

1.4 x 10" 10 

20, 21 

11. CH 3 + H0 2 — *CH 2 OH + OH 

2.6 x 10 ^ 

18 

12. CH, + HO, — CH . + O, 
3 2 4 2 

2.0 x 10" 12 exp-400/RT 

18 

13. CH , + CH,— CH. + CH, 
3 3 4 2 

4.0 x 10 ^ (Forward Only) 

18, 29, Text 

14. CH 3 + CH 2 0 — CH 4 + CHO 

2.0 x 10" 14 T 0,5 exp-600/RT 

18 

15. CH 3 + CHO— CH 4 + CO 

5.0 x 10' 13 T 0,5 

18 

16. CH 2 + 0 2 — CO + H 2 0 

4.0 x 10" 11 exp-2, 000/RT 

30, Text 

17. CH 2 OH + M — CH 2 O + H + M 

2.0 x 10" 9 exp-29, 000/RT 

20 

is. ch 2 oh + o 2 — ch 2 o + ho 2 

1.7 x 10" 11 

20 

19. CH 2 0 + M— CHO + H + M 

1.4 x 10' 7 exp-81, 000/RT 

19 

20. CH z O + OH — CHO + H 2 0 

3.8 x 10" 15 T 1 ' 4 

18, 22, Text 

21. CH 2 0 + H— CHO + H 2 

1.7 x 10“ 12 T°* 5 exp-3, 300/RT 

23, 33, Text 
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Table B-l (Continued) 


22. CH z O + O — CHO + OH 

3.0 x 10' 11 exp-3, 100/RT 

24 

23. CH z O + HO z — CHO + H 2 0 2 

1.7 x 10‘ 12 exp-7,900/RT 

18, 20 

24. CHO + M— CO + H + M 

4.0 x 10 exp-29, 600/RT 

19, 20 

25. CHO + O, —CO + HO, 

3.0 x 10* 13 T 0,5 

20,25,26, 

2 2 
26. CHO + OH — CO + H 2 0 

2.0 x 10 -10 

Text 
18, 19 

27. CHO + H — CO + H 2 

5.0 x 10 _1 ° 

25 

28. CHO + O —CO + OH 

2.0 x 10" 10 

27 

29. CHO + H0 2 — CO + H 2 0 2 

2.0 x 10 -I ° exp-3, 000/RT 

18 

30. CHO + CHO — CH 2 0 + CO 

2.0 x 10" 10 exp-700/RT 

25, Text 

31. CO + 0 2 — -CO z + O 

3.7 x 10‘ 13 exp-60, 000/RT 

28 

32. CO + OH — C0 2 + H 

1.1 x 10' 13 exp(T/H00) 

28 

33. CO + H0 2 — C0 2 + OH 

2.2 x 10' 10 exp-23,000/RT 

28 

34. CO + O + M — C0 2 + M 

4.0 x 10' 33 exp-4, 300/RT 

9 


C. Silane Oxidation Mechanism 
(To be used with and CH^ mechanisms) 


Reaction (Reversible) 


3 « 1 
Rate Coe££icient, cm -particle - 

-1 

sec unit 


1. SiH. + M — *-SiH, + H, + M 

4 2 2 

2. SiH . + O, — SiH,0 + H + OH 

4 2 £- 

3. SiH . + OH — SiH, + H,0 

4 3 2 

4. SiH. + H —SiH, + H, 

4 3 2 

5. SiH. + O — SiH, + OH 

4 3 

6. SiH. + HO, — SiH, + H,0, 

4 Z 3 Z Z 


2.5 x 10 2 T' 6 2 exp-59, 600/RT 

7.0 x 10' 8 exp-34,000/RT 

1.4 x 10" 11 exp-100/RT 

2.4 x 10'** exp-2, 500/RT 

7.0 x 10“ 12 exp- 1, 600/RT 

1.0 x 10" 11 exp-2, 000/RT 
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Table B-l (Continued) 


7. SiH 3 + C> 2 — SiH z O + OH 

8. SiH 3 + SiH 3 — SiH 4 + SiH 2 

9. SiH 3 + OH — •* SiH 2 + H £ 0 

10. SiH 3 + H — *-SiH 2 + H 2 

11. SiH 3 + O — SiH 2 <D + H 

U. SiH 3 + H0 2 -^SiH 4 + 0 2 

13. SiH 3 + SLH 2 0— -SLH 4 + SiHO 

14. SiH 3 + SiHO — SiH 4 + SiO 

15. SiH 2 + 0 2 —♦SiHO + OH 

16. SiH z O + OH —".SiHO + H z O 

17. SiH 2 0 + H — ►SiHO + H 2 

18. SiH 2 0 + O— SiHO + OH 

19. SiH z O + H0 2 — ►SiHO + H 2 0 2 

20. SiHO + M — ►SiO + H + M 

21. SiHO + 0 2 —SiO + H0 2 

22. SiHO + OH —SiO + H z O 

23. SiHO + H —SiO + H 2 

24. SiHO + O —SiO + OH 

25. SiHO 4- H0 2 — SiO + H 2 O z 

26. SiHO + SiHO — SiH 2 0 + SiO 

27. SiO + 0 2 — ►Si0 2 + O 

28. SiO + OH— Si0 2 '+ H 

29. SiO + H0 2 — Si0 2 + OH 

30. SiO + O + M— Si0 2 + M 


1.0 X 10" 10 exp-9, 700/RT 

1, Text 

2.0 x 10" 10 

T ext 

1.4 x 10~ U 

Text 

2.4 x 10" 11 exp- 1, 100/RT 

Text 

1.4 x 10" 10 

Text 

2.0 x lO -12 exp-400/RT 

Text 

2.0 x 10" 14 T 0,5 exp-600/RT 

Text 

5.0 x 10" 13 T 0,5 

Text 

4.0 x 10" 11 exp-2, 000/RT 

Text 

4.0 x 10 13 T * 

Text 

2.0 x 10" 12 T 0,5 exp-3, 300/RT 

Text 

3.0 x 10" 11 exp-3, 100/RT 

Text 

2.0 x 10" 12 exp-7, 900/RT 

Text 

1.0 x 10" 9 exp-32, 000/RT 

Text 

3.0 x 10" 13 T 0,5 

Text 

2.0 x 10" 10 

Text 

5.0 x 10' 10 

Text 

2.0 x 10“ 10 

Text 

2.0 x 10~ 19 exp-3, 000/RT 

Text 

2.0 x 10" 10 exp-700/RT 

Text 

1.0 x 10 13 exp-10000/RT 

Text 

1.1 x 10" 13 exp(T/ 1100) 

Text 

2.2 x 10" 10 exp-23, 000/RT 

Text 

4.0 x 10" 33 exp-4, 300/RT 

Text 
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Table B-l (Concluded) 


31. SiO z + SiO z -*-SL0 2 X + Si0 2 X 4.0 x 10' 8 T" 2 

32. SiH 2 0 + SiH 2 0 — SiH 2 OX 4.0 x 10" 8 T" 2 

+ SiH 2 OX 


Text 

Text 


D. Methane-Silane Oxidation Mechanisms Interaction Reactions 


Reaction (Reversible) 


1. SiH 4 + CH 3 — SiH 3 + CH 4 

2. SiH 3 + CH 3 —SiH 2 + CH 4 

3. SiH 3 + CHO — SiH 2 + CH z O 

4. SiH z O + CH 3 — SiHO + CH 4 

5. SiH z O + CHO— SiHO + CH z O 

6. SiHO + CH 3 — SiO + CH 4 

7. SiHO + CHO —SiO + CH 2 0 

8. SiO + CHO —SiHO + CO 


3 -1 

Rate Coefficient, cm -particle 
-1 

sec units 

Reference 

1.3 x 10' 12 exp-7, 000/RT 

40 

1.0 x 10' 12 exp-500/RT 

Text 

5.0 x 10' 13 T 0,5 

Text 

2.0 x 10“ 14 T 0,5 exp-600/RT 

Text 

2.0 x 10" 10 exp-700/RT 

Text 

5.0x10' 13 T°- 5 

Text 

2.0 x 10' 10 exp-700/RT 

Text 

2.0 x 10" 10 exp-700/RT 

Text 


Condensed Species; Si0 2 X, SiH 2 OX 

M Body Catalytic Weighting Factors; 

AR = 0.4, N 2 = 1.0, 0 2 = 1.1, H 2 = 2.0, H z O = 12.0 

H 2 0 2 = 16.0, CO = 1.0, CO z = 3.0, CH 4 = 1.5, 
CH z O = 6.0, CH 3 OH = 9.0, SiH 4 = 2.5, 

SIO = 16.0, Si0 2 = 16.0; 

All Free Radicals = 16.0 
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two-dimensional code and can readily accommodate mixing between streamlines 
and/or accept prescribed area flow channels. It is anticipated that these 
options will be exercised in subsequent efforts. 

Thermophysical input data were taken preferentially from the NBS JANAF 
Thermochemical Tables for species tabulated therein, or from the unpublished 
input data tapes supporting the computations of Ref. B-l, prepared by A.G. 
McLain of NASA- Langley Research Center for species (such as Sit^O, SiHO) 
not included in the JANAF tables. Exceptions to this are the standard heat 
of formation of SiH^ (46.6 kcal/mole), Sii^ (57.9 kcal/mole) and SiH 
(90.2 kcal/mole) which were obtained from the recent work of Doncaster and 
Walsh (Ref. B-4), Si02X - silica solids - thermo physical properties uti- 
lized are those corresponding to JANAF's "high cristobalite. " SiH^OX is a 
solid postulated to account for evidence of hydrogen bonding and alternate 
solids forms observed by early investigators of silane oxidation (e.g., 
Emeleus and Steward, Ref. B-5), with thermophysical inputs which have been 
estimated based on those for SiC^* SK^X and Si^O. 

In view of: (1) the present rather large ignorance factor of the indi- 

vidual species M-body weighting factors for use particularly in the silane 
reactions, and (2) the need to reduce the reaction mechanism's complexity 
for computational economy, a simplified approach was taken in this study and 
a single set of generic M-body weighting factors was utilized. These are 
given in Table B-l, relative to ^-which was assigned unit weight. They 
are derived from an interesting approximately linear correlation between log 
k and the boiling temperature of M, recently reviewed in the survey kinetics 
text by Kondratiev and Nikitin (Ref. B-6, p. 110). Here an upper asymptotic 
limit of 16 has been inferred from the scant data base. In any event, with 
^ the main constituent in the flows of interest, with rates entered for 
N 2 , results should not be overly sensitive to a reasonably judicious 
choice of generic M-body weighting factors, while the resulting simplifica- 
tion in the computations is appreciable. 
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A. Hydrogen Oxidation Mechanism 


For the most part, the hydrogen mechanism is that used in Ref. B-l, 
taken from Slack and Grillo (Ref. B-7), minus N0 x reactions. Rate coeffi- 
cient expressions for reactions 2 through 5 are, however, from the recent 
critical review by Cohen and Westberg (Ref. B-8)*; similarly, the rate 
coefficient for reaction 6 is from the recent NBS data compilation by 
Hampson (Ref. B-9). Additionally, reactions 15 through 17 have been added 
to the mechanism based on the finding of the sensitivity study by Dougherty 
and Rabitz (Ref. B-10) that they are required; rate coefficients for these 
reactions are those recommended by these authors. Similarly, the O+O+M 
reaction has been deleted from the mechanism, again on the basis of the 
sensitivity study by Dougherty and Rabitz. 

B. Methane Oxidation Mechanism 


As is evident from the references cited in Table B-lb, several methane 
kinetics sources were utilized, preferentially those from critical reviews 
by Westbrook and Dryer (Ref. B-18), Gardiner and Olson (Ref. B-19), and 
Tsuboi and Hashimoto (Ref. B-20). The methane mechanism is however less 
well settled than that for hydrogen and major discrepancies still exist 
between the assessments of the various authors. It was consequently found 
necessary to judiciously attempt to resolve some of these discrepancies and 
also to provide some additional data updates. These only will be briefly 
discussed below. 


The pre-exponential factor for reaction 2 has been modified (well within the 

-14 

uncertainty bound cited by Cohen and Westberg) from 7.5 x 10 to 1.0 x 
-13 3 -1 -1 

10 cm - particlec -sec , for purposes of obtaining closer correspondence 
between the computations made with the mechanism of Ref. B-l and this work. 
(At 1000K, the ratio of rates for this important reaction is 0.75 for the 
present set of rates versus that of Ref. B-l; forcing closer agreement would 
be less consistent with the error bounds cited by Cohen and Westberg in 
Ref. B-8. 
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Reaction 1. Tabayashi and Bauer's (Ref. B-ll) result for M-body dis- 
sociation of methane by Argon has been converted to the corresponding value 
with N£ as the M-body, utilizing the generic M-body weighting factors 
employed in this study. Agreement with Westbrook and Dryer's (Ref. B-18) 
expression at 2000K is excellent. This modification has been made also for 
other reactions involving M bodies, where appropriate. 

Reaction 2, The rate coefficient for the CH^ + C >2 reaction is 
poorly known. Huffington et al, (Ref. B-12) cite the "general belief" that 
O 2 attack is dominant below 2000K. Here, the pre-exponential factor has 
been estimated at one-twentieth that of the corresponding H 2 +O 2 reaction 
and the activation energy has been set equal to that for the H 2 +O 2 re- 
action. The resulting rate coefficient is about one-half that of the CH^ 

+ M rate (i.e.. Reaction 1) at 2000K and rapidly dominates it at lower tem- 
peratures due to the appreciably lower activation energy. A posteriori, the 
reasonableness of the results obtained in this study suggest that the pres- 
ent estimate is acceptable in the absence of definitive data. 

Reaction 5. The rate expression developed by Roth and Just (Ref. B-16) 
has been modified by this author to better accommodate the prior data and 
the more recent higher temperature data of Felder and Fontijn (Ref. B-17). 

Reaction 8. As observed in the review by Gardiner and Olson (Ref. 

B-19) there is considerable uncertainty in the temperature behavior of the 
rate of the CH^ + O 2 branching reactions, particularly at lower tempera- 
tures. Here, the approach has been taken to use the rate coefficient for 
Reaction 7 to represent the high temperature branch of the CH^ + O 2 re- 
action and to develop the rate coefficient for the low temperature branch of 
the identical reaction as Reaction 8, using the temperature dependence sug- 
gested by Bhaskaran et al. (Ref. B-31), fitting the pre-exponential factor 
to extensive classical ignition delay data correlated by Asaba et al. (Ref. 
B-32). We performed iterative calculations of ignition delay defined as the 
time required to reach 5 percent of the final equilibrium temperature rise 
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at constant pressure, as in Ref. B-l. It was found that the resulting delay 
times were almost inversely proportional to the sum of the rates for Reac- 
tions 7 and 8, with Reaction 8 dominating at low temperature as expected. 
Excellent agreement was obtained between the shock tube ignition delay times 
and the delay times so calculated, using the entire methane (and hydrogen) 
reaction mechanism as given in Table B-l, at 1 atm and initial temperatures 
varying from 700 to 1500K. However, a more proper computation for compar- 
ison to shock tube data would utilize a constant density (rather than pres- 
sure) post-shock condition. To do this with the ALFA code, a few code 
revisions need to be made which unfortunately have not yet been made. We 
hope to do this In the near future in an extension of this effort and repeat 
the aforementioned iterative computations to further refine the mechanism. 
For this reason, these computed results will not be presented here. It is 
important to note however that the changes resulting from the computational 
refinements are not anticipated to be more than by a factor of 2 at most, 
based on prior experience and other data analysis. Thus while the changes 
are thought to be important and necessary for proper kinetic interpretation 
of shock tube data, they are not of overriding importance and do not negate 
the basic goodness-of-fit of results obtained with the present analysis. 

Reaction 9. Tsuboi and Hashimoto's (Ref. B-20) rate coefficient for 
the CHg + OH reaction has been employed, with the reaction written to 
yield Cl^OH + H, as favored by Gardiner and Olson (Ref. B-19). 

Reaction 10, Washida's (Ref. B-21) room temperature result for this 
reaction has been used, with the zero activation energy indicated by Tsuboi 
and Hashimoto (Ref. B-20). 

Reaction 13. The radical sink reaction between CH^ + CH^ is known 
to be important in the methane mechanism (see, e.g., Ref. B-18). The major 
product is actually thought to be 02 ^, rat ^ er than CH^ + CHj as 
written in Table B-l to avoid complexities arising from the introduction of 
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C 2 hydrocarbons into the mechanism. For this reason we allow this global- 
ized reaction only to proceed in the forward direction in the computations, 
using the rate coefficient measured by Baulch and Duxbury (Ref. B-29) and 
thus dodging unnecessary complexities in the mechanism. 

Reaction 16. The room temperature result of Laufer (Ref. B-30) has 
been used for the C^ + reaction together with an estimated activa- 
tion energy of 2 kcal/mole. To eliminate spurious chain branching from the 
effective radical sink reaction CH^ + CH 3 , here written as leading to 
CH^ + CH 2 , the products of the CH 2 + O 2 reaction have been written 
as CO + H 2 O, rather than CHO + OH. Thus the combined result of Reactions 
13 and 16 is, effectively, CH 3 + CH 3 + 0 2 -*■ CH 4 + CO + H 2 0, i.e., 
a chain termination step. 

Reaction 20. Low temperature recent results by Stief et al. (Ref. 

B-22) have been pooled with higher temperature results cited by Westbrook 
and Dryer (Ref. B-18) to obtain the rate coefficient expression used in this 
study. 

Reaction 21. Klemm's (Ref. B-23) low temperature expression for the 
rate coefficient of the C^O + H reaction has been here rewritten to ob- 
tain agreement with higher temperature results of Dean, Johnson and Steiner 
(Ref. B-33). 


Reaction 25. The rate coefficient given in Table B-l for the CHO + 

O 2 reaction is the result of reconciliation by the author of low tempera- 
ture data by Reilly et al. (Ref. B-25) and Shibuya et al. (Ref. B-26) with 
the high temperature result of Tsuboi and Hashimoto (Ref. B-20). 

Reaction 30, The rate coefficient given in Table B-l for the CHO + CHO 
reaction is the result of assigning a pre-exponential factor equal to that 
for the reaction of CHO with 0, OH and HO 2 (Reactions 26, 27, and 28) and 
equating an Arrhenius rate coefficient format to the room temperature result 
of Reilly et al. (Ref. B-25), thereby obtaining a quite reasonably activa- 
tion energy of 700 cal/mole. 
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C. Silane Oxidation Mechanism 


Beyond the most elementary reactions of SiH^ itself with M (thermal 
dissociation), OH, H and 0, virtually none of the detailed steps of the 
silane combustion mechanism have been firmly established, much less studied 
to the point where measurements of individual rate coefficients are avail- 
able. Consequently this mechanism is largely in an exploratory status. The 
present formulation has utilized the early mechanism derived by Beach, 
Mackley, Rogers and Chinitz (Ref. B-l) as a starting point, with the advan- 
tage of the first experimental shock tube measurements of silane-hydrogen- 
oxygen-nitrogen delay times, by McLain (Ref. B-34). Only a very brief 
discussion of the basis of the mechanism is permitted here. 

Reaction 1. In Ref. B-l the products of the thermal dissociation reac- 
tion SiH^ + M were written as SiH^ + H, based largely on an analysis of 
the results of Strater (Ref. B-35). Strater's work however was performed in 
a reactor packed with silica chips and it seems likely that the results 
therein were largely dominated by heterogeneous processes. This conclusion 
is strongly reinforced by the demonstration by Baliga and Ghandi (Ref. B-36) 
that silane decomposition rates under conditions similar to Strater's are 
strongly dependent on substrate temperature and obey the rate theory of 
heterogeneous reactions at a solid surface. Moreover, the existing homoge- 
neous phase thermal decomposition data of Newman et al. (Ref. B-37) and also 
of Purnell and Walsh (Ref. B-38) is quite convincing in the finding that 
H£ rather than H is the hydrogenic species resulting from the primary 
thermal dissociation step. For these reasons this author has felt compelled 
to accept the results of Refs. 37 and 38, writing the thermal decomposition 
step as SiH^ + M — ^SiH 2 + H 2 + M, This of necessity required addi- 
tional extensive revision to the mechanism of Ref. B-l as discussed below. 

Reaction 2. Silane is well-known to react spontaneously - and also 
somewhat erratically - with oxygen at temperatures near ambient over a wide 
range of equivalence ratios. This can only be possible - regardless of an 
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apparent role of impurities of either gas phase or heterogeneous origin as 
chain starters - if a direct reaction between 0^ and SiH^, as postulated 
in Table B-l, occurs. The rate coefficient listed is the result of fitting 
early ALFA computations of ignition delay times - defined as discussed in B 
above - to the recent shock tube measurements of ignition delay of McLain 
(Ref, B-34). Better agreement with the shock tube data and some minor re- 
vision of the rate coefficient is anticipated when ALFA can model a constant 
density flow regime to more accurately describe conditions behind the re- 
flected shock used in NASA’s ignition experiments. 

Reactions 3, 4, and 5. The rate coefficients for the elementary reac- 
tions of silane with OH, H and 0 are those recommended in the study and 
review of Atkinson and Pitts (Ref. B-39) and Arthur and Bell (Ref. B-40), 
respectively. 

Rate coefficients for the remaining reactions are all estimates, some 
of which have had the benefit of refinement possible because of the avail- 
ability of McLain's recent data (Ref. B-34). 

Reaction 6. The Arrhenius parameters for the SiH^ + HO 2 abstrac- 
tion reaction have been estimated from those for the corresponding reaction 
of silane with OH, H, and 0 (i.e.. Reactions 3, 4, and 5). 

Reaction 7. The rate coefficient for the SiH^ + O 2 reaction, 
written as in Ref. B-l, is that estimated by Chinifz (Ref. B-l), with the 
activation energy modified to 9700 cal/mole based on observation by 
Vasilyeva et al. (Ref. B-41) of this activation energy for the upper igni- 
tion limit ratio of silane to O 2 at atmospheric pressure. This assignment 
must of course be regarded as tentative until such time as a deeper analysis 
of the available silane explosion limit data is attempted. Such analysis 
unfortunately was beyond the scope of this limited study. 

Reaction 8. Heavy reliance has been placed on the logic employed in 
developing the methane mechanisms discussed above to write the products of 
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the radical sink SiH^ + SiH-j reaction as SiH^ + SiH 2 , with a 
temperature-independent rate coefficient, quite in analogy to the CH^ + 

CH-j radical sink reaction. The rate coefficient has been increased to 2.0 
x 10 ^ from the 4,0 x 10 ^ of the CH^ + CH^ reaction partly on an 
intuitive basis and partly as a result of mechanistic refinement required to 
better match the McLain data. 

Reactions 9 and 10. The rate coefficients for the SiH^ + OH and 
SiH^ -I- H abstraction reactions have been estimated as follows: Preexpo- 

nential factors equal to those for the corresponding SiH^ reactions have 
been assigned; the activation energy for the OH reaction is estimated to be 
less than that for the OH reaction wih SiH^, i.e., essentially zero; the 
activation energy for the H atom reaction has been estimated at 1100 cal/ 
mole based on the corresponding rate coefficient for the H abstraction re- 
action with disilane reported by Arthur and Bell (Ref. B-40) and the convic- 
tion that the SiH^ rate is appreciably faster than the SiH^ rate. 

Raction 11. The SiH^ + 0 reaction has been written as the complete 
analogue of the CH^ + 0 reaction, with identical rate coefficient, in the 
absence of any other input. 

Raction 12. Similarly, the SiHg + H0 2 + SiH^ + 0 2 reaction 
has been written as the silane analogue of the corresponding CHj + 

H0 2 CH^ + 0 2 reaction, with identical rate coefficient. Note that 
the analogue of the CH^ + H0 2 CH 2 0H + OH reaction has been omitted 
from the mechanism to avoid having to treat the largely unknown S1H 2 0H 
species. 

Reaction 13. Again, the CH^ + CH 2 0 -*■ CH^ + CH0 reaction has 
served as the postulated analogue for the postulated SiH^ + SiH 2 0 -*■ 

SiH^ + SiHO reaction, with identical rate coefficient. 
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Reaction 14, Here the analogue has been the CH^ + CHO CH^ + CO 
reaction, with presumed essentially equal rate coefficient for the postu- 
lated SiH 3 + Si HO -*• SiH^ + SiO reaction. 

Reaction 15. Again, the logic employed in the methane system for the 
CH 2 + O 2 reaction has been invoked for the Si^ + O 2 analogue, with 
an equal rate coefficient assigned. Note that as in the methane mechanism, 
the only Si^ reaction entered is the dominant one with O 2 . 

Reactions 16 through 30. With one exception all of these postulated 
reactions have been written as proceeding just as their analogues in the 
methane mechanism, with essentially equal or slightly faster rate coeffi- 
cients. The exception is the SiO + 0 2 -* Si02 + 0 reaction which has 
been written with a rate coefficient which has Relatively little impact on 
ignition delay times at higher temperatures but without which no ignition 
can occur at temperatures near ambient even with quite appreciable sensiti- 
zation by trace species. With the mechanism as written, auto-ignition of 
silane-air mixtures at temperatures near ambient has been tentatively ident- 
ified as being due to sensitization by such trace species, as from contami- 
nants and/or the products of heterogeneous reactions. For example, ozone in 
amounts often present at sea level (i.e.’, 1 to 10 parts per billion ) is pre- 
dicted to induce ignition of stoichiometric silane-air mixtures at 1 atm and 
300 K in from 2 to 7 sec with the mechanism as given in Table B-l. With no 
ozone or similarly reactive free radical-generating agent - or without 
Reaction 30 - no perceptible reaction occurs. 

Reactions 31 and 32. Condensation of solids results in appreciable 
(i.e. , several hundreds of degrees Kelvin) temperature rise in silane com- 
bustion, the magnitude of course depending on the concentration of silane. 
Inasmuch as such large temperature changes drastically alter rates for indi- 
vidual reactions - especially those with high activation energy - accurate 
kinetic modeling of silane ignition requires inclusion of solids nucleation. 
Here, a very simple, expedient approach to modeling solids condensation 
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kinetics has been taken: Rather than becoming embroiled in the exact de- 
tails of the very complex nucleation processes - largely unknown for sila- 
ceous species at this time in any event - a simple ternary condensation 

process, e.g., Si0 9 + Si0 9 + M ■+ Si0 9 X + Si0 9 X + M was invoked, with a 

z z z -30 -1 

rather ordinary ternary rate coefficient of about 6 x 10 T Further, 

since the ALFA code cannot handle such a reaction type (without special 

modification which seemed unwarranted) the reaction was rewritten as a 

simple binary rate process with rate coefficient k x [M] , i.e., with [M] 

evaluated at atmospheric pressure, the rate process is, e.g., Si0 9 + SiO, 

fg _2 ^ 

-> Si02X + S102X with a rate coefficient approximated by 4.0 x 10 T , 
quite adequate for the relatively small pressure excursion range 0.5 to 2.0 
atm contemplated for application of the mechanism developed for this study. 

D . Me thane-Silane Oxidation Mechanisms Interaction Reactions 


This subset extends the silane mechanism developed above for use with 
the hydrogen and methane mechanisms simultaneously, as required for modeling 
silane-assisted ignition of methane. Note that the reactions have all been 
written in the exothermic direction. 

Reaction 1. Arthur and Bell (Ref. B-40) reviewed results for the SiH^ + 
CHj -»■ SiHg + CH^ reaction, and their recommended value has been utilized 
in the mechanism developed. 

Reaction 2. The rate coefficient for the presumed CH^ + SiH^ CH^ + 
SiH 2 reaction has been estimated from the not disparate values for the CH^ + 
CH 2 O CH^ + CHO and CH^ + HO 2 CH ^ + O 2 reactions, in the absence of 

any other input data. 

Reactions 3, 4, 6, and 7. These processes fyave been written exotherm- 
ically, with rate coefficients equal to the equivalent encounters in the 
methane system, i.e., the CH^ + CHO encounter for SiH^ + CHO and the 
CH 2 O + CH^ encounter for Si^O + CHg, etc. 
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Reactions 5 and 8. The Sii^O + CHO and S10 + CHO processes have aLso 
been written exothermically, each with rate coefficients equal to that for 
the CHO + CHO encounter, in the absence of any other input data. 

E. Status of Mechanism 


In synopsis, it is clear that where experimental inputs have not been 
available, very heavy use of insights and data inputs developed for the 
methane combustion mechanism have been imbedded in this work into the devel- 
opment of the silane combustion mechanism, much as was the principle of the 
approach taken by Chinitz (Ref. B-l) in developing the earlier silane mech- 
anism. Here the methodology has been carried one step further, giving the 
present mechanism benefit of a large body of literature pertaining to meth- 
ane combustion where required to close gaps in the relatively meager silane 
data base. The extent to which this adaptation of the methane formalism 
accurately describes silane results is still a moot point. As the silane 
data base expands additional mechanistic refinement can surely be antici- 
pated. 


Particular kinetic issues of concern at the moment include: (1) the 

need for Independent experimental confirmation of such key rates as that 
postulated between silane and 0 2 (i.e., Reaction 2, Table B-lc); (2) the 
need for a detailed exploration of the sensitivities of predicted combustor 
heat release profiles to variations in uncertain rate coefficients and minor 
reaction pathways, leading to greater understanding of the SiH^/H^CH^ 
system and, also, to a mechanism reduced to its essentials; (3) the extent 
of reaction, if any, between SiO and the species H 2 0 and C0 2 . As pres- 
ently constituted (with the SiO + OH reaction the primary step leading to 
Si0 2 ) in fuel-rich regimes with excess H 2 in the flow, the present mech- 
anism leads to an undershoot in flame temperature potentially of several 
hundred degrees Kelvin, depending on mixture composition. This has been 
traced to incomplete (i.e,, nonequilibrium) combustion of SiH^ to SiO due 
to removal of uncombined 0 2 by rapid combustion of H 2 in excess, 
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effectively stopping the SiO oxidation chain. Is this a true nonequilib- 
rium, cool flame phenomenon? Does SiO react with I^O, or CO 2 , to form 
Si02? Russian work indicates this possibility - as well as the possi- 
bility that perhaps a solid SiO phase needs to be considered: "Solid SiO 

reacts slowly with atmospheric O 2 at room temperature. Complete oxidation 
is not achieved at 500 C. It is pyrophoric in a finely divided state and 
burns to Si02 in air, producing a flame. Water vapor reacts noticeably 
with SiO even at 500 C, producing Si02 and I^. Silicon monoxide is 
slowly oxidized in a CO 2 atmosphere at 400 C...” (Berezhnoi, (Ref. B-42, 
p. 102), and; (4) finally, definitive determination of the identity and con- 
centration history of reaction intermediates and end-products - gaseous and 
solid - is required to assure the understanding of silanp combustion chem- 
ical physics. 

RESULTS 


Digests of results computed in this study have been forwarded under 
separate cover at periodic intervals to Dr. R.C. Rogers, Hypersonic Propul- 
sion Branch, Langley Research Center. A detailed summary of these extensive 
results is not presented in the present report due to time and space limita- 
tions. We will however summarize the high points and principal results of 
our effort. 

A. Silane-Hydrogen Combustion 


Figure B-l shows results obtained for the recent shock tube experiments 
of McLain (Ref. B-34). The uppermost curve shows results computed by McLain 
usin the reaction mechanism of Ref. B-l, employing the NASA code described 
in Ref. B-43 run in a constant density flow option. T ignition t * ie 
ignition delay time based on extrapolation of a T versus time arithmetic 
coordinate plot at its steepest slope back to the temperature of the initial 
post-shock gas mixture. Also shown on the figure is McLain's experimental 
data and the computed results obtaind in the present study with the mechanism 
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Species 


Mole Fraction 



1000/T.K 

Fig. B-l - Ignition Times for 20/80 Percent Si lane/ Hydrogen 
P 0 = 1.25 atm (All Calculations) 
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of Table B-l and the ALFA code (Ref. B-3) run in a constant pressure option, 

using t . , as defined above, and also, t T he latter delay time 
° Ignition 5%AT 

corresponds to the time required to reach 5 percent of the fully equili- 
brated temperature rise, as in Ref. B-l. As already observed elsewhere in 
this text, a more correct calculational procedure for analysis of shock tube 
ignition delay measurements involves computation of T ig n ition a con ~ 
stant density flow option, as by McLain - which option was not available in 
the particular version of the ALFA code employed in this study. However, 
based on experience, computations of x ig n -Ltj_ on at constant density fall 
somewhere between the computations of and T ig n i t j; on at constant pres- 

sure. Thus the curves computed using the ALFA code and the present reaction 
mechanism, which straddle the McLain data to within at least a factor of two 
are also expected to straddle computed values of T jg n iti 0n at constant 
density, using the present reaction mechanism. Agreement with the present 
data is thus to within at least a factor of two. 


B. Silane-Methane Combustion 


Parametric computations of methane combustion in air flows with silane 
as combustion aid have been performed with the refined mechanism. These 
computations extend over Initial temperatures of from 600 to 1200 K, at con- 
stant pressures from 0.5 to 2.0 atm, and overall stoichiometries from 0.25 
to 1.0, with from 2 to 20 vol.% silane in the methane fuel. Consistent with 
early, unpublished engine tests with hydrocarbons at GASL, we find that use 
of silane as an ignition aid materially reduces required ignition delays and 
combustion times for methane in air. For example, at 1 atm constant pres- 
sure and an initial temperature of 1000 K, the computed value of x^g^^^ 

is reduced from 0.70 sec for stoichiometric CH. - air by three orders of 

-4 

magnitude - to 5.3 'x 10 sec - with 10 vol.% SiH^ in the fuel; the cor- 
responding time requirement to reach 95% of the equilibrium combustion tem- 

-4 

perature rise is reduced from 0.73 to 9.2 x 10 sec, i.e., also by nearly 

three orders of magnitude. For the range of parameters explored, 

and x 95 *£ ^ are found to vary essentially as the first power of reciprocal 
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(constant) pressure. The T versus time plots show interesting structure 
with the mechanism as presently constituted. Typically there is a rela- 
tively steep, accelerating initial temperature rise which gradually becomes 
less steep as reactions proceed (i.e., showing some cool flame character- 
istics in early stages) and then gradually increases again, with the slope 
characteristically being steepest at T ig n i t i on ' At this point the cause 
of this very interesting and unexpected behavior has not been ascertained. 

CONCLUSIONS 


Detailed kinetics mechanisms for silane-hydrogen and methane oxidation 
have been combined and adapted for use in computation of the characteristics 
of supersonic air-methane-silane/I^ ignition and heat release. The re- 
sulting formulation results in a workable combustion mechanism which (with 
the ALFA code) typically requires only nominal computational time - say 5 to 
20 CRU for a representative 95% heat release run, with the NASA Langley Re- 
search Center CDC 6600 computer. 

After initial mechanistic refinement interfacing with very recently ob- 
tained NASA-Langley shock tube ignition data, predictions made with the re- 
fined mechanism agree to at least within a factor of 2 with these early 
shock tube measurements of ignition delay over the range 800 to 1000 K for 
silane-H2~02/N2 mixtures. Additionally, a sensitivity to ignition pro- 
motion in silane-air mixtures at temperatures near ambient by trace species 
arising from contaminants and/or heterogeneous processes has been demon- 
strated with the present mechanism, again in keeping with experimental ob- 
servation. 

Initial favorable results for use of silane as an ignition aid in 
supersonic methane combustors have been obtained in a parametric examina- 
tion of theoretical ignition and reaction times in silane-CH^-air mixtures. 
Reductions in ignition and reaction times by factors on the order of 3 
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orders of magnitude is predicted, depending on particulars. Use of silane 
as an ignition aid thus appears to offer considerable promise for augmented 
combustion of hydrocarbon fuels, as well as for hydrogen. 


RECOMMENDED NASA SILANE PROGRAM DIRECTIONS 


Based on the favorable results of this study and our interaction with 
and understanding of the status of the NASA silane hypersonic combustion 
program, we offer the following recommendations for the further implementa- 
tion and study of this very promising advanced combustion technique: 


• Measure silane-hydrogen/oxygen-nitrogen ignition delays in 
shock tube experiments over a range of initial tempera- 
tures, pressures and compositions, utilizing existing com- 
putational kinetics codes to support analysis of the ex- 
periments. Analysis of combustion products and measurement 
of combustion temperature, particularly in fuel-rich mix- 
tures with excess hydrogen, is required to assess the role 
of nonequilibrium combustion. (The latter work could per- 
haps be accomplished using a flat flame burner to supple- 
ment the shock tube.) 

• Extend the measurements and analysis to silane-hydrocarbon 
combustion, first to methane, then to propane, and ulti- 
mately to highly practical flight candidates such as JP- 
type aviation fuels. 

• Aggressively pursue engine tests using silane as an ig- 
niter, as well as a pilot, for each of the main fuels of 
interest (i.e., H2, CH4, C3H8, JP, etc.). Couple these 
engine tests with kinetic analysis of the combustor flow 
fields, including recirculation effects known to be of 
paramount importance to successful flameholding. (A three- 
dimensional code is required, necessitating use of combus- 
tion mechanisms reduced to their essentials only.) 

• Explore the use of silane itself, or of higher silanes, as 
the primary fuel in a hypersonic propulsion system. 
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Appendix C 

GLOBALIZATION OF HYDROCARBON KINETICS 
FOR HYPERSONIC COMBUSTION COMPUTATIONS 



Appendix C 


I 


f 


SUMMARY 


A generalized semi-global model of hydrocarbon flame front kinetics 
suitable for use with the GIM code has been developed for analysis of hyper- 
sonic combustion. Real fuels of most interest are liquid hydrocarbon 
aircraft-type fuels, and the model development has accordingly specifically 
addressed these fuels. 
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C.l BACKGROUND 


Modeling the details of combustion of real hydrocarbon fuels including 
the intricacies of complex pre-ignition kinetics remains remote for two 
primary reasons: 

1. The mechanistic details of the combustion of higher hydro- 
carbons are still largely unresolved. 

2. Even if the detailed kinetics were unambiguously known, the 
magnitude of the computational problem would still swamp 
existing computational facilities. 

To overcome this, partially empirical combustion models for hydro- 
carbons ("global, semi-global and quasi-global” models, recently reviewed by 
Chinitz in Ref. C-l) are under development - all of which assume the hydro- 
carbon to decompose in a single, fictitious global reaction to relatively 

simple products, viz. CO, 1^0 and/or C^H^ species such as CH^ or 
C 2 H 4* 


Global models utilize irreversible global rate expressions derived from 
curve fits to data for all of the combustion products. For example, the 
global model of Dryer et al. (Ref. C-2) for propane combustion utilizes the 
following four global reactions to describe the temporal behavior of the 
species C^Hg, ° 2 > ^0, CO 2 and I^O: 

C 3 H 8 .+ 1.5 C 2 H 4 + H 2 
C 2 H 4 + 0 2 + 2 CO + 2 H 2 
CO + 0.5 0 2 -> C0 2 
H 2 + 0.5 0 2 -*■ H 2 0 

None of the individual steps in such models bears any semblance to the 
fundamental free radical-dominated kinetic steps painstakingly documented in 
studies of actual combustion mechanisms undertaken over the course of this 
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century. This results in rather obvious shortcomings and great uncertainty 
in extrapolation beyond the narrow data range fit by the computationally 
straightforward, but restrictive, expressions employed. 

Semi-global models relax this major shortcoming by allowing CO/I^/ 

H2O intermediate reaction products to follow their individual, detailed 
fundamental free radical-dominated kinetic p^ths while utilizing a global 
expression for consumption of the hydrocarbon to intermediate species, e.g., 
as in the global, irreversible finite rate step proposed by Harsha et al. 
(Ref. C-3) and recommended by Chinitz (Ref. C-l): 


C 

n 


H m + T °2 


2. H 2 + n CO 


Quasi-global models allow an additional degree of complexity by incorp- 
orating intermediate yields of C^-C 2 hydrocarbon decomposition products. 
Individual detailed fundamental kinetic paths of these species are then 
followed along with reactions of the C0/H 2 /H 2 0 products of the initial 
globalized hydrocarbon pyrolysis reaction (e.g., C H 2n+2 + M * ( n /2) 

C 2 H, + H, + M) . 

At this time, none of the models is clearly best and all have both 
shortcomings and controversial features outside of the limited regimes for 
which they have been developed. In particular, existing global models do 
not give acceptable results in flows dominated by long ignition delay times. 

After review of the problem we concluded that the best available 
reasonably computationally efficient yet data-based compromise lay in devel- 
opment of a generalized semi-global model of flame front and post-flame 
front kinetics with heavy reliance on experimental measurements of ignition 
delay for the specific hydrocarbon considered. 
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C. 2 TECHNICAL APPROACH 


Experimental measurements of auto-ignition delays (kinetically-limited) 
are usually available for fuels of interest over the range of parameters 
likely to be encountered. If not available, experimental facilities exist 
which can generate the required data. Such data are to be utilized heavily 
in the construction of a relatively simple but realistic model of the com- 
bustion kinetics of (generic) hydrocarbon fuels, as follows. Coupled with 
fundamental CO/^/^O rate processes the resulting model is adequate for 
finite rate analysis of hypersonic combustion by the GIM or other complex 
flowfield modeling codes. 


The temporal development of combustion in a given parcel of fuel-air 
mixture is conceptually divided into three regions: 

Region 1 - Physical Processes 


In this region pre-ignition processes are dominated by physical 
processes which include droplet formation, heating, vaporization, diffusion, 
mixing and temperature equilibration with the air stream. Typically, this 
region is less important than Region 2, below. Region 1 is characterized by 
a required physical processes residence time, T . 


Region 2 - Auto- Ignition Kinetics 


In this region pre-ignition processes are dominated by the highly com- 
plex chemical processes which result in the partial decomposition of high 
molecular weight hydrocarbon species and the formation of critical concen- 
trations of intermediate free-radical species, with only a small change in 
temperature and in the concentration of major species (i.e., fuel and C^). 
Region 2 is characterized by a required autQ-ignition residence time, T^. 
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Region 3 - Flame Front and Post-Flame Front Kinetics 


In Region 3, hydrocarbon molecules are rapidly consumed by oxidation, 
resulting in a hot flame front with accompanying combustion of l^/CO/ 

C ^-02 intermediates (C^-C^ species are not incorporated in the present 
model) followed by post-flame front free radical recombination processes. 

Figure C-l summarizes these regions and Fig. C-2 summarizes results of 
several workers for auto-lgnltlon delays of aircraft-type fuel sprays in 
air, with auto-ignition delays for Hj superimposed. 

METHODOLOGY 

Region 1. Pre-ignition physical processes (i.e., spray and droplet 
formation, vaporization, mixing, temperature equilibration, etc.) are 
assumed to be separately analyzed. Normally the physical time scale will be 
a multiple - a - of the time required for the initial liquid jets - moving 
at velocity v^ - to traverse a characteristic cross channel flow dimension 
X, where X will as often as not be a spray bar half spacing, or its equiv- 
alent. Thus, 


T - (X/v.) a (C.l) 

P L 

where a is a number presumed of order unity obtained in ancilary computa- 
tions of the time requirement for these physical events preceding ignition. 

All reaction rates are set equal to zero in Region 1. 

Region 2 . Auto-ignition delay data like those shown in Fig. C-2 are 
assumed to be available for a specific fuel-air combustion problem. Func- 
tionally is input as 


N 


l AI 


a ai T 
P b . f(<j>) 


e ai 

exp Rf" sec 


(C.2) 
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Fig. C-2 - Autoignition of Liquid Hydrocarbon Fuel Sprays 

in Air (Largely ffom Ref. Cr4, as cited in Ref. C-9) 


r . I - 22 » i° — ej 

AI, Hydrocarbons qP 

8.0 x 10 9600 

t ai,h 2 " p exp T“ 


sec (Ref. 


(Boundary) 

C-12) 


C-7 




where 


* Pre-exponential constant of curve fit to auto-ignition 
data 

T = Temperature in Region 2, K 

N ■ Temperature exponent (default * 0) 

P = Pressure in Region 2, atm 

b • Pressure exponent (default ■ 1) 

f ( 4> ) * Function of fuel-air stoichiometry ratio, (default 

value of f(<j>) * 1) 

_ (Fuel; Air molar ratio) 

^ ” (Fuel : Air molar ratio at stoichiometry) 

E A i ■ Auto-ignition delay apparent activation energy, cal/mole 
R = Gas constant 


All reaction rates are set equal to zero in Region 2 also, correspond- 
ing to the pre-ignition condition of hardly noticeable change in T, P and 
major species concentration profiles evident in plots such as those of Fig. 
A-l and A-2. The onset of ignition at the end of Region 2 - or rather the 
beginning of Region 3 - is modeled by inputting a discrete, characteristic 
post-ignition value for the mole fraction of chain branching H atom free 
radical - - and initiating the globalized kinetics as discussed below. 

Region 3 


The kinetics are switched on ip Region 3, with the mole fraction of H 
atoms initialized at a level corresponding to post-ignition, pre-flame con- 
ditions, i.e., 
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Input value 


H, 


H, 


-4 -3 

A number of order 1 x 10 to 1 x 10 
(see e.g., Figs. A-l and A-2). 


(C.3) 


Note that the characteristic rise time for the H atom population to 

increase from X^ to a level typical of flame fronts, i.e., a number of 
_2 “I 

order 10 is a fraction of the autoignition delay for H», i.e., x AT u . 

t A1 y tl2 

From Fig. C-2 it is evident that ^ is appreciably shorter than 

T AI, Hydrocarbons* 111118 the H atom rise tlme from \ to flame front 
levels is only a relatively small fraction of the total auto-ignition delay 

for hydrocarbon combustion. Results are accordingly relatively insensitive 

-4 

to the choice of X . For hydrocarbon combustion, 1 x 10 is a good 

H I _3 

choice for X . For H_ combustion, 1 x 10 is however more appro- 

H I 1 _4 _3 

priate as the corresponding H atom rise time from 10 to 10 mole 


fraction can be a significant fraction of T 


AI, H 2 * 


A rate coefficient - - for depletion of hydrocarbon at the en- 

trance to Region 3 is input such that the time constant for hydrocarbon 
oxidation is equal to HC /10. Thus, the hydrocarbon concentration is 
decreased by finite rate reaction to 1/e of its initial value in a time 
equal to 10 percent of T ^ ^ - i.e., the flame front thickness is set at 
about 10 percent of the kinetic ignition length. Thus, letting 


CH+ y 0, n CO + y H 
n m l z l l 


(C.4) 


be the irreversible finite-rate oxidation reaction for the hydrocarbon, with 
rate coefficient k^ defined by 

d [C H ] 

" — at" m - t°n V t° 2 l < C - 4a > 
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and 


d [0 2 ] 


dt 

- ' T “hc l C n »„) i° 2 l 

(C.4b) 

d [CO] 
dt 

■ n “hC > C n »„! I° 2 ] 

(C.4c) 

d [h 2 ] 

' T^IC t C „V [° 2 1 

(C.4d) 

dt 

the rate coefficient is given by 


k HC ^ °2 

J _ 1 10 

o t hc t ai,hc 

(C.5) 

Thus 




10 


[( V T AI,HC 
o ’ 


- ■ 

10 • pt • <«) 

[° 2 J A I» 
o 

(C.6) 

For air 



[0 2 ] - (P/T) 

x 1.54 x 10 21 particle cm -3 

(C.7) 


0 


with P in atm and T in K, neglecting fuel dilution. Thus 


"hc 


6.5 x 10' 21 \ b-1 _~ E AI 

P . f (<J>) . exp-g^- 


A A! T 


N-l 


(C.8) 


3 , -1 -1 
cm . particle . sec 
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To fix this concept, consider the aircraft fuel ignition delay data of 
Fig. C-2. While each individual data set can be represented by a specific 
expression, if this is desired, all of the data are bounded by the expression 


1.22 x 10" 6 8400 

T AI,HC = P ex P sec 


as indicated in the caption to Fig. C-2. Employing this as representative 
of a generic aircraft fuel hydrocarbon, substitution in Eq. (C.8) yields the 
expression 

k H( , = 5.3 x 1(T 15 T exp-16, 700/RT (C.9) 

Default values for N, b and f(<J>) have been used in the above expression 
in the absence of more specific inputs, which may or may not be available 
for a particular case. 

The reaction set is completed by adding those elementary reactions 
necessary and sufficient to model l^/CO afterburning combustion effi- 
ciently as in the basic plume chemistry models utilized in Lockheed's 
analyses of SSME/SRB/ Trident, etc., advanced propulsion systems. The com- 
plete semi-global hydrocarbon kinetics model is shown in Table C-l. Ele- 
mentary rate coefficients for the reversible reactions 2 through 9 are all 
from Table B-l and have previously been discussed. 


Table C-l - SEMI-i 
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Appendix D 

PROGRESSIVE ASSEMBLY OF THE GIM CODE ELEMENTS 



Appendix D 


D.l INTRODUCTION 


This appendix describes the GIM code difference analogs from a slightly 
different perspective than that developed in the previous documentation. 

This development is based on new weight functions which allow the integra- 
tion of the analogs in a closed form which is amenable to vectorization and 
repetative calculation. As a result, the analogs can now be computed with 
great efficiency when needed (i.e., "progressively") rather than occupying 
machine storage. This saving in storage is particularly important in the 
solution of large three-dimensional problems. 

It is assumed here, as in all GIM development, that all functions of 
the local variables can be expressed in a discretizing relation as inter- 
polated sums of point functions: 

f(n) - s 1 (n) f i5 i-i,...,2 d 

where d is the problem dimensionality. Further, the shape functions, , 
are assumed here to be multilinear Lagrange interpolants. 

The following summation conventions hold, except where noted: 

1. Repeated indices are summed. 

2. Lower case Latin indices run over the 2^ corner points. 

3. Greek indices run over the d dimensions. 

4. Capital Latin indices run over all nodes in the discrete 
domain. 
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D.2 TWO-DIMENSIONAL ANALOGS 


The two-dimensional divergence law equation is 


3U + 3E + 3F 
3t 3 x 3 y 


0 


This is modeled by the discrete difference equation on the element e (Fig. 
D-la): 


where 


and 


! 6 = " ° 


'ij 


'ij 


'ij 


= f do. f 1 do, w. s . 
J i y 2 i j 


O O 


3 3 3(s 1t y) 

" / dn i / dn 2 "i 


o. o 


J- J- 3(x,S.) 

/ dn, ydn 2 W t 


o o 


si - (i - m) (i - n 2 ) 

S 3 - 01 02 


S 2 - m (1 - 02) 
S 4 - (1 - 01 ) 02 
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I 



1 


a. Element Domain 


I 

e = 2 

1 - o 

e - 1 

e » 3 

► ' < 

e * 4 

>. ■ ■■ ' ■■ < 


b. Assembled Computational Domain 


Fig. D-l - Two-Dinjensional Rectilinear Domains 


I 
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The choice of the following quadratic weight functions: 


w x - (9 - 36 n x + 30 nj) (9 - 36 n 2 + 30 n 2 ) 
W 2 - (3 - 24 + 30 nj) (9 - 36 n 2 + 30 n 2 ) 
W 3 - (3 - 24 n x + 30 nj) (3 - 24 n 2 + 30 n 2 ) 
W A - (9 - 36 n 1 + 30 nj) (3 - 24 n 2 + 30 n 2 ) 


yield closed form analogs of the differential expressions. These can be 
written in the following form: 

e 


ij 


6(x 1 , y ± ) 

0 


i “ j 
i “ i 


B ij E j 


'<S(E i , y i ) 
6 ( 11 !. 0 2 ) 


c « ' 




n 2 ) 

where the finite difference Jacobian determinant is defined as 


6(11!* ri2^ 


0 0 

C. A, A S-i 
Xv X 1 v 1 


with no sum on i. 


e is the two dimensional Levi-Civita symbol and A 

AV V 


f^ is the two point finite difference in the direction evaluated at 
point i in element e. 
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The element equations can be assembled over the four elements which 
make up the typical 9-node computational domain by performing the weighted 
Boolean sum 


^ a e l& " ° 
e«l 


where 

4 

E*. ‘ 1 

e*l 

This results in an assembled field equation for node N 


a nm u m 


+ B. 


NM 




+ c 


NM ‘M 


0 


where 




4 

Ea 


NM — : ’ e 
e»l 


6 <V y N> 

5 (n ^ 9 H 2 ^ 




NM 


e«l 


^ ^ E N* y N^ 

6(ri 1 > n 2 > 


C NM F M " a e 


g ( x N* V 

SCrii > n 2 ) 


Since Aj^ is diagonal, the resulting difference scheme is an explicit 
finite difference analog to the divergence law. 
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D.3 THREE-DIMENSIONAL ANALOGS 


The approach to the three-dimensional analogs is directly analogous to 
the 2-D derivation. The three-dimensional divergence law equation 


9U 3F 3G 

3t 3x 3y 8z 


is modeled by the element equation 


<1 u j + b ij e j + c h f j + D u 


where e is the eight-node element (Fig. D-2a) and 


ij 


/ d 3 n w A Sj 


9(x,y,z) 

3(n 1 .n 2 »n 3 ) 


ij 


- f d 3 n w A 


3(Sj,y,z) 

3 (n 3 »fl2»^3^ 


'ij 


f d 3 n w i 


3 (x»Sj ,z) 

3 ("n »n 2 * r > 3 ) 


'ij 


/ d3 " w i 


3 (x>y>Sj) 

3 (ri *ti 2 3 ^ 


(i - n,) (l - n 9 ) (l - n,) 


(l - t )^ (i - n 2 ) n 3 


O-i (1 — Ho) (i “ Oo) 


n l (1 ■ V n 3 


» 3 ■ n x n 2 (i “ n 3 ) 


S 7 " ’ll ^2 ^3 


(i - np n 2 (i - n 3 ) 


5 8 " (1 " ’V n 2 n 3 
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Choosing the following weight functions: 

w x = 64 (4 - 30^ + 60nJ - 35nJ) (4 - 3on 2 + 60n 2 - 35n 2 ) 

. (4 - 30n 3 + 60n 3 - 35n 3 ) 

W 2 ■* 64 (- 1 + 15^ - 45nJ + 35n 3 ) (4 - 30n 2 + 60n 2 - 35 ti 2 ) 

. (4 - 30n 3 + 60n 3 - 35n 3 ) 

w 3 = 64 (- i + isn 1 - 45nJ + 35n 2 ) (- i + i5n 2 - 45n 2 + 35n 2 ) 

. (4 - 3on 3 + 60n 3 - 35n 3 ) 

W 4 - 64 (4 - 30ri 1 + 60nJ - 35nJ) (- 1 + i5n 2 - 45n 2 + 35n|) 

. (4 - 30n 3 + 60n 2 - 35n 3 ) 

W 5 - 64 (4 - 301-^ + 6onJ - 35nJ) (4 - 30n 2 + 60n 2 - 35 ti 2 ) 

. (- 1 + 15n 3 - 45H^ + 35n^) 

W 6 - 64 (-1 + I5n 1 - 45nJ + 35nJ) (4 - 3on 2 + 6on* - 35n 2 ) 

. (- l + i5n 3 - 45n 3 + 35n 3 ) 

W 7 - 65 (- i + I5n 1 - 45nJ + 35nJ) (- i + I5n 2 - 45 n^ + 35n^> 

. (- 1 + I5n 3 - 45n 3 + 35n 3 ) 

w g - 64 (4 - 30n 1 + 6onJ - 35 nJ) (- 1 + i5n 2 - 45n 2 + 35n 3 ) 

. (- 1 + i5n 3 - 45n 2 + 35n 3 ) 
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element analogs result which are similar to the two-dimensional form, i.e., 


A 6 . 


= < 


6(n 1 *n2*n 3 ) 


i * j 


B" E. = 
ij 3 


<S(E i> y 1 ,z i ) 

S(n 1 ,n 2 ,n 3 ) 


ct. f. = 

1 3 3 


6(n 1 »n 2 * r i 3 ) 


DT. G. = 
ij 3 


6 (x^y^G^ 
6 (g ! »r| 2 ,r l 3 ) 


Here the three-dimensional finite difference Jacobian determinant is 


6(^,8^^) 

6(n 1 »n 2 *n 3 >. 


e e e 

e, A. f, A g, A h. 
Ayv A i y & i v i 


(no sum on i). 

The three-dimensional assembly consists of a weighted Boolean sum over 
the eight elements which comprise the classic 27 node finite difference cell 
(Fig. D-2b). 

The assembled equation becomes 

A X7W U + B X7W EL, + C xlVf F w + D XTU G m - 0 
NM M NM ^1 NM M NM M 


where 




8 

^ a A® 
e-1 


E«. 


e °NM 


»n 2 »n 3 ) 
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NM 


e=l 


<(E H’ y W r B ) 

^ > ^2 > ^ 3 ) 


8 

C F m = 2-f a 

NM M e 

e=l 


^^ x n >f n ,2 n^ 
6 > ^2 > 


d nm g m 


= Zoi, 


e=l 


5< x N >y N >G N ) 

6 (n^ » r l2 ,r t3 ) 


8 

it a e = 1 

e=l 


D . 4 SUMMARY 


This appendix has presented an alternative development of the GIM 
finite difference analogs, A different choice of weight functions provides 
closed-form discrete difference expressions of the integral analogs. These 
expressions are amenable to vectorization and repetative calculation. By 
calculating the analogs progressively during the solution procedure much 
storage can be saved. 


There are a number of other, less obvious, advantages to this develop- 
ment : 

• Fewer nodal connections are necessary (5 in 2-D; 7 in 3-D). 

• No significant differences have been seen in comparisons 
with earlier GIM code versions. 

• The MATRIX module is no longer necessary 

• The GEOM module performance is significantly improved by 
eliminating the element matrix integration 


D-10 



• The INTEG module is now Independent of the GEOM module. 

Any grid generation program can now be used so long as its 
output is suitably formatted. 

• The finite difference scheme can now be changed without 
multiple GEOM runs. 
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